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Abstract 

We describe a subdivision algorithm for isolating the complex roots 
of a polynomial F € C[a;]. Given an oracle that provides approximations 
of each of the coefficients of F to any absolute error bound and given an 
arbitrary square B in the complex plane containing only simple roots of F, 
our algorithm returns disjoint isolating disks for the roots of F in B. 

Our complexity analysis bounds the absolute error to which the co¬ 
efficients of F have to be provided, the total number of iterations, and 
the overall bit complexity. It further shows that the complexity of our 
algorithm is controlled by the geometry of the roots in a near neighborhood 
of the input square B, namely, the number of roots, their absolute values 
and pairwise distances. The number of subdivision steps is near-optimal. 
For the benchmark problem, namely, to isolate all the roots of a polynomial 
of degree n with integer coefficients of bit size less than r, our algorithm 
needs 0{n^ + n^T) bit operations, which is comparable to the record bound 
of Pan (2002). It is the first time that such a bound has been achieved using 
subdivision methods, and independent of divide-and-conquer techniques 
such as Schonhage’s splitting circle technique. 

Our algorithm uses the quadtree construction of Weyl (1924) with two 
key ingredients: using Pellet’s Theorem (1881) combined with Graeffe 
iteration, we derive a "soft-test" to count the number of roots in a disk. 
Using Schroder’s modified Newton operator combined with bisection, in 
a form inspired by the quadratic interval method from Abbot (2006), 
we achieve quadratic convergence towards root clusters. Relative to the 
divide-conquer algorithms, our algorithm is quite simple with the potential 
of being practical. This paper is self-contained: we provide pseudo-code 
for all subroutines used by our algorithm. 
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1 Introduction 


The computation of the roots of a univariate polynomial is one of the best 
studied problems in the areas of computer algebra and numerical analysis, 
nevertheless there are still a number of novel algorithms presented each year; 
see [26, 24, 25, 27, 33] for an extensive overview. One reason for this development 
is undoubtedly the great importance of the problem, which results from the 
fact that solutions for many problems from mathematics, engineering, computer 
science, or the natural sciences make critical use of univariate root solving. 
Another reason for the steady research is that, despite the huge existing literature, 
there is still a large discrepancy between methods that are considered to be 
efficient in practice and those that achieve good theoretical bounds. For instance, 
for computing all complex roots of a polynomial, practitioners typically use 
Aberth’s, Weierstrass-Durand-Kerner’s and QR algorithms. These iterative 
methods are relatively simple as, in each step, we only need to evaluate the given 
polynomial (and its derivative) at certain points. They have been integrated 
in popular packages such as MPSolve [5, 6] or eigensolve [16], regardless of 
the fact that their excellent empirical behavior has not been entirely verified in 
theory. In contrast, there exist algorithms [15, 28, 32] that achieve near-optimal 
bounds with respect to asymptotic complexity; however, implementations of 
these methods do not exist. The main reason for this situation is that these 
algorithms are quite involved and that they use a series of asymptotically fast 
subroutines (see [32, p. 702]). In most cases, this rules out a self-contained 
presentation, which makes it difficult to access such methods, not only for 
practitioners but also for researchers working in the same area. In addition, for 
an efficient implementation, it would be necessary to incorporate a sophisticated 
precision management and many implementation tricks. Even then, there might 
still be a considerable overhead due to the extensive use of asymptotically fast 
subroutines, which does not show up in the asymptotic complexity bounds but 
is critical for input sizes that can be handled on modern computers. 

In this paper, we aim to resolve the above described discrepancy by presenting 
a subdivision algorithm for complex root isolation, which we denote by CIsolate. 
For our method, we mainly combine simple and well-known techniques such as 
the classical quad-tree construction by Weyl [57], Pellet’s Theorem [39], Graeffe 
iteration [4, 19], and Schroder’s modified Newton operator [49] . In addition, 
we derive bounds on its theoretical worst-case complexity matching the best 
bounds currently known for this problem; see Section 1.1 for more details. Hence, 
we hope that our contribution will finally bring together theory and practice 
in the area of complex root finding. In this context, it is remarkable that, for 
the complexity results, we do not require any asymptotically fast subroutines 
except the classical fast algorithms for polynomial multiplication and Taylor shift 
computation. Our presentation is self contained and we provide pseudo-code 
for all subroutines. Compared to existing asymptotically fast algorithms, our 
method is relatively simple and has the potential of being practical. 

In theory, the currently best algorithm for complex root finding goes back to 
Schonhage’s splitting circle method [48], which has been considerably refined by 
Pan [32] and others [21, 31] . In [32], Pan gives an algorithm for approximate 
polynomial factorization with near-optimal arithmetic and bit complexity.^ From 

^Pan considers a similar model of computation, where it is assumed that the coefficients of 
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an approximate factorization, one can derive isolating disks for all complex roots. 
A corresponding algorithm for complex root isolation, which uses Pan’s method 
as a subroutine, has been presented and analyzed in [28]. Its cost can be 
expressed in terms of (accessible) parameters that directly depend on the input 
such as the degree of F and the size of its coefficients, but also in terms of 
(hidden) geometric parameters such as the pairwise distances between the roots. 
A special case, namely the so-called (complex) benchmark problem of isolating 
all complex roots of a polynomial F with integer coefficients of bit size at most 
r, has attracted a lot of interest in the literature. Using Pan’s method [15, 28], 
the latter problem can be solved with 0{n^T) operations^, which constitutes the 
current record bound for this problem.^ So far, there exists no other method for 
complex root isolation that achieves a comparable bound. For the real benchmark 
problem, that is the isolation of the real roots of a polynomial of degree n with 
integer coefficients of bit size at most r, recent work [45] describes a practical 
subdivision algorithm based on the Descartes method and Newton Iteration 
with bit complexity 0{n^ + n^r). An implementation of this method [22] is 
competitive with the fastest existing implementations ]41] for real root isolation, 
and it shows superior performance for hard instances, where roots appear in 
clusters. Our contribution is in the same line with [45], that is, both methods 
combine a subdivision approach, a simple predicate to test for roots, and Newton 
iteration to speed up convergence. The main difference is that we treat the more 
general problem of isolating all complex roots, whereas the algorithm from ]45] 
can only be used to compute the real roots, due to the use of Descartes’ Rule of 
Signs to test for roots. 

We further remark that, in comparison to global approaches such as MP- 
SOLVE ]5, 6], which compute all complex roots in parallel, our algorithm can also 
be used for a local search for only the roots contained in some given square. In 
this case, the number of iterations as well as the cost of the algorithm adapt to ge¬ 
ometric parameters that only depend on the roots located in some neighborhood 
of the given square. 


the input polynomial are complex numbers that can be accessed to an arbitrary precision. Then, 
for a polynomial F with roots 21 ,..., 2:71 contained in the unit disk and an integer L > n log n, 
Pan’s algorithm computes approximations Zi of Zi with IjP’—lcf(F)-]~[ (LiU-iOlli <2-^-||F||i 
using only O(nlogL) arithmetic operations with a precision of 0{L). For a lower bound on 
the bit complexity of the approximate polynomial factorization, Pan considers a polynomial 
whose coefficients must be approximated with a precision of Q(L) as, otherwise, the above 
inequality is not fulfilled. This shows that already the cost for reading sufficiently good 
approximations of the input polynomial is comparable to the cost for running the entire 
algorithm. Hence, near-optimality of his algorithm follows. In the considered computational 
model. Pan’s algorithm also performs near-optimal with respect to the Boolean complexity 
of the problem of approximating all roots. However, we remark that this does not imply 
near-optimality of his method for the benchmark problem of isolating the complex roots of an 
integer polynomial. Namely, Pan’s argument for the lower bound is based on a lower bound 
on the precision to which the coefficients have to be approximated. In the case of integer 
polynomials, the coefficients are given exactly, hence the cost for reading an arbitrary good 
approximation of the polynomial never exceeds the cost for reading the integer coefficients. 

^With 0(-}, we indicate that poly-logarithmic factors are omitted, i.e., for a function p, we 
denote with 0{p) the set of functions in 0{plog^ p), where c is a constant. 

^So far, the bound 0{n^T) can only be achieved by running Pan’s factorization algorithm 
with an L of size r2(n(r - 1 - logn)), which means that ©(n^r) bit operations are needed for any 
input polynomial; see [15, Theorem 3.1] for details. The adaptive algorithm from [28] needs 
0{n^ -h n^r) bit operations, however its cost crucially depends on the hardness of the input 
polynomial (e.g., the separations of its roots), hence the actual cost is typically much lower. 
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1.1 Overview of the Algorithm and Main Results 

We consider a polynomial 

n ^ 

F{x) = ^ a, X* S C[x], with n>2 and - < |a„| < 1. (1) 

i=0 

Notice that, after multiplication with a suitable power of two, we can always 
ensure that the above requirement on the leading coefficient is fulfilled, without 
changing the roots of the given polynomial. It is assumed that the coefficients 
of F are given by means of a coefficient oracle. That is, for an arbitrary L, the 
oracle provides a dyadic approximation hi of each coefficient Ui that coincides 
with Oi to L bits after the binary point. We call an approximation F obtained 
in this way an (absolute) L-bit approximation of F and assume that the cost for 
asking the oracle for an L-bit approximation of F is the cost of reading such an 
approximation;"^ see Section 2 for more details. Let us denote by zi to the 
roots of F, where each root occurs as often as determined by its multiplicity. 
Now, given a closed, axis-aligned square B in the complex plane, our goal is to 
compute isolating disks for all roots of F contained in B. Since we can only ask 
for approximations of the coefficients, we need to further require that B contains 
only simple roots of F as, otherwise, a multiple root of multiplicity k cannot be 
distinguished from a cluster of k nearby roots, and thus the problem becomes 
ill-posed. If the latter requirement is fulfilled, then our algorithm CIsolate 
computes isolating disks for all roots contained in B.^ However, it may also 
return isolating disks for some of the roots contained in 2B, the square centered 
at B and of twice the size as B. Our approach is based on Weyl’s quad tree 
construction, that is, we recursively subdivide B into smaller sub-squares and 
discard squares for which we can show that they do not contain a root of F. 
The remaining squares are clustered into maximal connected components, which 
are tested for being isolating for a single root. 

As exclusion and inclusion predicate, we propose a test based on Pellet’s 
theorem and Graeffe iteration. We briefly outline our approach and refer to 
Section 3 for more details. Let A ;= A(to, r) C C be the disk centered at m with 
radius r, and define A • A(to, r) := A(m, A • r) for arbitrary A € K’*'. According 
to Pellet’s theorem [39], the number of roots contained in A equals k if the 
absolute value of the fc-th coefficient of La(x) := F{m + rx) dominates the sum 
of the absolute values of all other coefficients. For k = 0 and fc = 1, it has been 
known [46, 59] that Pellet’s theorem applies if the smaller disk • A contains 
k roots and the larger disk • A contains no further root, where ei and 62 are 
suitable positive constants. In the paper at hand, we derive constants ei and 
62 such that the latter result stays true for all fc. As a consequence, using only 
O(loglogn) Graeffe iteration for iteratively squaring the roots of Fa, we can 
replace the factors and by the constants pi := « 0.94 and p 2 '■= |. 

^Notice that we only require approximations of the coefficients, hence our method also 
applies to polynomials with algebraic, or even transcendental coefficients. In any case, the 
given bounds for the cost of isolating the roots of such a polynomial do not encounter the cost 
for computing sufficiently good L-bit approximations of the coefficients. Depending on the 
type of the coefficients, this cost might be considerably larger than the cost for just reading 
such approximations. 

®If the requirement is not fulfilled, our algorithm does not terminate. However, using an 
additional stopping criteria, it can be used to compute arbitrarily good approximations of all 
(multiple) roots; see the remark at the end of Section 4.2 for more details. 
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More precisely, we derive a test that allows us to exactly count the number of 
roots contained in a disk A, provided that the disks p 2 ■ A and pi ■ A contain the 
same number of roots. If the latter requirement is not fulfilled, the test might 
return the value —1, in which case we have no information on the number of 
roots in A. Since, in general, the latter test requires exact arithmetic and since 
we can only ask for approximations of the coefficients of F, there might be cases, 
where we either cannot decide the outcome of our test or where an unnecessarily 
high precision is needed. Based on the idea of so-called soft-predicates [61], we 
formulate a variant of the above test, which we denote by T*, that uses only 
approximate arithmetic and runs with a precision demand that is directly related 
to the maximal absolute value that F takes on the disk A. 

In the subdivision process, we inscribe each square in a corresponding disk 
and run the T»-test on this disk. Squares, for which the test T, yields 0, do not 
contain a root and can thus be discarded. The remaining squares are clustered 
into maximal connected components, which we also inscribe in corresponding 
disks. If the T*-test yields 1 for such a disk, we discard the cluster and store the 
disk as an isolating disk. Otherwise, we keep on subdividing each square into 
four equally sized sub-squares and proceed. This approach on its own already 
yields a reasonably efficient algorithm, however, only linear convergence against 
the roots can be achieved. As a consequence, there might exist long paths in the 
subdivision tree with no branching (there are at most n — 1 branching nodes). 
For instance, when considering the benchmark problem, there exist polynomials 
(e.g. so called Mignotte polynomials having two roots with a very small distance 
to each other) for which the length of such a sequence is lower bounded by fl(nT). 
We show how to traverse such sequences in a much faster manner, that is reducing 
their length to 0(log(nr)) in the worst-case, via a regula falsi method, which 
combines Newton iteration and square quartering. Our approach is inspired by 
the so-called quadratic interval refinement (QIR for short) method proposed 
by Abbott [1]. It combines the secant method and interval bisection in order 
to further refine an interval that is already known to be isolating for a root. 

In [42, 43, 45[, the QIR approach has been considerably refined by replacing the 
secant method by Newton iteration (i.e. Schroder’s modified Newton operator 
for multiple roots). Compared to Abbott’s original variant, this yields a method 
with quadratic convergence against clusters of roots during the isolation process. 
Our approach is similar to the one from [45[, however, we use the T»-test instead 
of Descartes’ Rule of Signs, which only applies to real intervals. Furthermore, 
the approach from [45[ uses fast approximate multipoint evaluation [21, 23[ in 
order to determine subdivision points whose distance to the roots of F is not 
too small. This is needed to avoid an unnecessarily large precision when using 
Descartes’ Rule of Signs. For our algorithm CIsolate, there is no need for (fast) 
approximate multipoint evaluation. We now state our first main theoretical 
result, which shows that our algorithm performs near-optimal with respect to 
the number of produced squares: 

Theorem. Let F he polynomial as in (1) and suppose that F is square-free. For 
isolating all complex roots of F, the algorithm CISOLATE produces a number of 
squares bounded by 

6 (n ■ log(n) • log (n • Ff • log(CTp^))) , 

where we define log(a:) := max(l, log |a;|) for arbitrary x € €.,Tp := lbg(max[h;^ \zi\) 
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the logarithmic root bound and ap := min(j \zi — Zj\ the separation of F. 

For the benchmark problem, the above bound simplifies to 0{n log(n) log(nT)). 
When running our algorithm on an arbitrary axis-aligned square B, we obtain 
refined bounds showing that our algorithm is also adaptive with respect to the 
number of roots contained in some neighborhood of B as well as with respect to 
their geometric location. Namely, suppose that the enlarged square 2B contains 
only simple roots of F, then we may replace n, Tp, and ap in the bound in 
the above theorem by the number of roots contained in the enlarged square 2B, 
the logarithm of the width of B, and the minimal separation of the roots of F 
contained in 2B, respectively; see also Theorem 6. 

Finally, we give bounds on the the bit complexity of our approach as well as 
on the precision to which the coefficients of F have to be provided: 

Theorem. Let F be a polynomial as in (1) and suppose that F is square-free. 
For isolating all complex roots of F, the algorithm CIsolate uses a number of 
bit operations bounded by 

O n- {Tp-\-n- lbg(2j) -k lbg(crF(2i)"^) -k lbg(i^'(2i)"^))) = 

0{n{n^ -I- nlog(MeaF) + lbg(Disc^^))), 

where we define Tp := [log ||F||ool; o'pizi) := \zi — Zj\ the separation of Zi, 

Meai? := |a„| max(l, \zi\) the Mahler Measure, and DisCi? the discriminant 

of F. As input, the algorithm requires an L-bit approximation of F with 

L = d + n ■ Ibg(zi) -f lbg(crF(^i)"^) + log(F'( 2 j)"^))^ 

= 0{n^ nlbg(MeaF) + lc)g(Disc[^^)). 

Again, we also give refined complexity bounds for the problem of isolating 
all roots of F contained in some square B, which show that the cost and the 
precision demand of our algorithm adapt to the hardness of the roots contained 
in a close neighborhood of the square. For the benchmark problem, the above 
bound simplifies to 0{n^ + n^r). It is interesting that our bounds on the bit 
complexity for isolating all complex roots as achieved by CIsolate exactly match 
the corresponding bounds for the complex root isolation algorithm from [28], 
which uses Pan’s method for approximate polynomial factorization. 

1.2 Related Work 

As already mentioned at the beginning, there exists a huge literature on com¬ 
puting the roots of a univariate polynomial. This makes it simply impossible 
to give a comprehensive overview without going beyond the scope of a research 
paper, hence we suggest the interested reader to consult some of the excellent 
surveys [26, 24, 25, 27, 33[. Here, we mainly focus on a comparison of our method 
with other existing subdivision methods for real and complex root finding. 

For real root computation, subdivision algorithms have become extremely 
popular due to their simplicity, ease of implementation, and practical efficiency. 
They have found their way into the most popular computer algebra systems, 
where they constitute the default routine for real root computation. Prominent 
examples of subdivision methods are the Descartes method [9, 13, 14, 41, 42, 
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45, 44, 51], the Bolzano method® [3, 7, 46], the Sturm method ]11, 12], and the 
continued fraction method ]2, 50, 53, 54]. From a high-level point of view, all of 
the above mentioned methods essentially follow the same approach: Starting from 
a given interval Iq, they recursively subdivide Iq to search for the roots contained 
in Iq. Intervals that are shown to contain no root are discarded, and intervals 
that are shown to be isolating for a simple root are returned. The two main 
differences between these algorithms are the choice of the exclusion predicate 
and the way how the intervals are subdivided. For the real benchmark problem of 
isolating all real roots of a polynomial of degree n with integer coefficients of bit 
size r or less, most of the above methods need 0{nT) subdivision steps and their 
worst-case bit complexity is bounded by O(n'^r^). The bound on the number of 
subdivision steps stems from the fact that the product of the separation of all 
roots is lower bounded by and that only linear convergence to the roots 

is achieved. By considering special polynomials (e.g., Mignotte polynomials) that 
have roots with separation one can further show that the bound 0{nT) 

is even tight up to logarithmic factors; see ]8, 14]. When using exact arithmetic, 
the cost for each subdivision step is bounded by 0(n®T) bit operations, which 
is due to the fact that n arithmetic operations with a precision of 0(n^T) are 
performed. In [44, 45], it has been shown for the Descartes method that it 
suffices to work with a precision of size 0{nT) in order to isolate all real roots, a 
fact that has already been empirically verified in ]41] . This yields a worst-case 
bit complexity bound of size 0(n®T^) for a modified Descartes method, which 
uses approximate instead of exact arithmetic. For a corresponding modified 
variant of the Bolzano method [3], a similar argument yields the same bound. 
Recent work ]42, 45, 51] combines the Descartes method and Newton iteration, 
which yields algorithms with quadratic convergence in almost all iterations. They 
use only 0(nlog(nT)) subdivision steps, which is near optimal. The methods 
from ]42, 51] work for integer polynomials only and each computation is carried 
out with exact arithmetic. An amortized analysis of their cost yields the bound 
0(n®T) for the bit complexity. ]45] introduces an algorithm that improves upon 
the methods from [42, 51] in two points. First, it can be used to compute the 
real roots of a polynomial with arbitrary real coefficients. Second, due to the 
use of approximate arithmetic, its precision demand is considerably smaller. For 
the real benchmark problem, it achieves the bit complexity bound 0(n® -I- n^r). 
More precisely, it needs 0(nlog(nT)) iterations, and, in each iteration, 0(n) 
arithmetic operations are carried out with an average precision of size 0{n + t). 
This essentially matches the bounds achieved by our algorithm CIsolate for 
complex root isolation. CIsolate shares common elements with the method 
from [45], however we had to develop novel tools to accommodate the fact 
that our search area is now the entire complex plane and not the real axis. In 
particular, we replaced Descartes’ Rule of Signs, which serves as the test for real 
roots in ]45], by our novel test T* for counting the number of complex roots in 
a disk. 

For computing the complex roots, there also exist a series of subdivision 
methods (e.g. ]10, 29, 30, 34, 38, 40, 46, 58, 60]); however, only a few algorithms 

®The Bolzano method is based on Pellet’s theorem (with k = 0). It is used to test an 
interval I for roots of the input polynomial F and its derivative F'. I contains no root if 
Pellet’s theorem applies to F. If it applies to F', the function F is monotone on I, and thus I 
is either isolating for a root or it contains no root depending on whether there is a sign change 
of F at the endpoints of I or not. 
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have been analyzed in a way that allows a direct comparison with our method. 
The earliest algorithm most relevant to our work is Weyl’s [57]. He proposed 
a subdivision based algorithm for computing a 2“^-relative approximation to 
all the roots of a polynomial, which is a slightly different problem then root 
isolation. The inclusion and exclusion tests are based on estimating the distance 
to a nearest root from the center of a box, or what are called proximity tests 
in the literature. The arithmetic complexity of the algorithm is 0{n^blogn), 
when not using asymptotically fast polynomial arithmetic. The problem with 
Weyl’s approach, indeed with any approach based on subdivision, is the linear 
convergence to the roots. The convergence factor was improved by Renegar 
[40] and Pan [34] by considering a combination of subdivision with Newton 
iteration. Renegar [40] uses the Schur-Cohn algorithm [18, Section 6.8] as 
an exclusion test (rather than the proximity tests of Weyl). In addition, he 
introduces a subroutine for approximating the winding number of a polynomial 
around the perimeter of some disk, and thus a method for counting the number 
of roots of the polynomial in a disk. Once the number k of roots in a disk is 
known, a fixed number (depending on the degree and the radius of the disk) of 
Newton steps are applied to the {k — l)-th derivative of the polynomial, which 
guarantees quadratic convergence to a cluster containing k roots. The arithmetic 
complexity of Renegar’s algorithm for the problem of approximating the roots is 
0{n^ log 6 + logn) without using asymptotically fast polynomial arithmetic. 
The improvement over Weyl’s result is basically due to the quadratic convergence 
obtained by the use of Newton iteration. 

Pan ]34] describes another modification of Weyl’s approach that has arithmetic 
cost 0((n^ logn) log(6n)), which is an improvement over Renegar’s algorithm 
since the dependence on the degree is a quadratic factor in n. The exclusion 
test is based on a combination of Turan’s proximity test [55] and Graeffe iter¬ 
ation. Note that the asymptotic complexity of these tests is 0{n), whereas a 
straightforward implementation of the Schur-Cohn test takes 0{n?) arithmetic 
operations; the difference in the cost of these exclusion tests is the reason behind 
the the improvement in the complexity estimate of Pan’s algorithm compared to 
Renegar’s. The algorithm in ]34] recursively interchanges Schroder’s iteration (a 
modification of Newton’s iteration to handle multiple roots) and Weyl’s subdivi¬ 
sion process. As in the case of Renegar, the former is needed to approximate a 
cluster of roots, and if that fails to happen, the subdivision is used to break up 
the set of roots into smaller subsets, and continue recursively. The transition 
between the iteration phase and the subdivision process is based on estimating 
the root radii [47, Section 14], and is perhaps more adaptive than Renegar’s 
approach. To estimate the number of roots inside a disk (which is needed to 
estimate the size of a cluster), Pan uses a combination of the winding number 
algorithm along with Graeffe iteration to ensure that there are no roots close 
to the boundary of the disk; as suggested by Pan, one can alternatively use 
the root radii algorithm without affecting the complexity significantly. The 
analysis of the algorithm has two steps. First, is to bound the number of boxes 
computed in the subdivision phase. This is done by considering the connected 
components of the boxes and bounding the number of boxes in each component 
in terms of the number of roots inside a slight scaling of the smallest disk 
containing the component; in our case, the bound on the number of boxes is 
obtained by mapping the components to appropriate roots (see Theorem 6); the 


resulting bound is comparable in both cases (see ([34, Prop. 8.3] in Pan and 
Theorem 5 below). The second step of the analysis shows that for certain well 
separated clusters Newton iteration gives us quadratic convergence to the cluster 
[34, Lem. 10.6]; an analogous result is also derived by Renegar [40, Cor. 4.5], 
and by us (Lemma 6). Some of the key differences between the approach in 
this paper and Pan’s [34] are the following: we use Pellet’s test combined with 
Graeffe iteration for both the exclusion test and detecting a cluster; we use a 
modification of the QIR method [1] for multiple roots, which is more adaptive 
in transitioning between the quadratic convergence and subdivision phases. In 
terms of the results derived, perhaps the most important difference is that we 
bound the bit complexity of our algorithm. In comparison, neither Renegar nor 
Pan analyze the precision demand or the Boolean complexity of their algorithms. 

Similar to our method, Yakoubsohn ]60] combines Weyl’s quad tree approach 
and a test for roots based on Pellet’s theorem. However, since only an exclusion 
predicate (based on Pellet’s theorem with fc = 0) is considered but no additional 
test to verify that a region is isolating, his method does not directly compute 
isolating regions but arbitrary good approximations of the complex roots. In ]46], 
we introduced a variant of Yakoubsohn’s method, denoted by Ceval, that 
computes isolating disks for the complex roots of an integer polynomial. There, 
an additional inclusion test (based on Pellet’s theorem with k = 1) has been 
used to show that a disk is isolating for a root. The methods from ]46, 60] only 
consider square-quartering, and thus nothing better than linear convergence can 
be achieved. For the benchmark problem, the algorithm from ]46] needs O^n^r) 
subdivision steps and its cost is bounded by 0{n‘^T^) bit operations. Yakoubsohn 
further mentions how to improve upon his method by combining the exclusion 
predicate with Graeffe iterations, which yields an improvement by a factor of size 
n with respect to the total number of produced squares. In [17], an extension 
of Pellet’s theorem for analytic functions has been considered and thoroughly 
analyzed. The authors also derive further criteria to detect clusters of roots of 
such functions, and to determine their multiplicities and diameters. This allows 
for the computation of suitable starting points for which Schroder’s modified 
Newton operator yields quadratic convergence to the cluster. In contrast, we 
follow the approach of combining Pellet’s theorem and Graeffe iteration to derive 
a simple test for detecting clusters of roots. However, we do neither compute 
the diameter of such a cluster nor do we consider any additional computations 
to check whether quadratic convergence to the cluster can be achieved. Instead, 
we rely on a trial-error approach that performs Schroder’s modified Newton 
operator by default and then checks for success. We show that this can be done 
in a certified manner such that quadratic convergence to clusters is guaranteed 
for all but only a small number of iterations, where our method falls back to 
bisection. Our approach works well with polynomials whose coefficients can only 
be approximated and we derive precise bounds on the precision demand in the 
worst-case. 

In our previous work [61], we provided the first complete algorithm for 
computing e-clusters of roots of analytic functions. Like the present work, it is a 
subdivision approach based on the T^-test of Pellet; but unlike this paper, it does 
not have quadratic convergence nor complexity analysis. In [61], we assumed 
that an analytic function is given when we also have interval evaluation of its 
derivatives of any desired order; this natural assumption is clearly satisfied by 
most common analytic functions. The algorithm from [61] does not compute 


9 


isolating disks but arbitrary small regions containing clusters of roots, hence 
being also applicable to functions with multiple roots and for which separation 
bounds are not known. 

1.3 Structure of the Paper and Reading Guide 

In Section 2, we summarize the most important definitions and notations, which 
we will use throughout the paper. We suggest the reader to print a copy of this 
section in order to quickly refer to the definitions. We introduce our novel test 
T* for counting the roots in a disk in Section 3. The reader who is willing to 
skip all details of this section and who wants to proceed directly with the main 
algorithm should only consider the summary given at the beginning of Section 3, 
where we give the main properties of the T*-test. The algorithm CIsolate 
is given in Section 4. Its analysis is split into two parts. In Section 5.1, we 
derive bounds on the number of iterations needed by our algorithms, whereas, 
in Section 5.2, we estimate its bit complexity. Some of the (rather technical) 
proofs are outsourced to an appendix, and we recommend to skip these proofs in 
a first reading of the paper. In Section 6, we summarize and hint to some future 
research. 


2 Definitions and a Root Bound 

Let f be a polynomial as defined in (1) with complex roots zi,..., Zn- We fix 
the following definitions and denotations: 

• As mentioned in the introduction, we assume the existence of an oracle that 
provides arbitrary good approximations of the coefficients. More precisely, 
for an arbitrary non-negative integer L, we may ask the oracle for dyadic 
approximations a-i = of the coefficients Oi such that rrii & X-\-i ■ X a <C 
are Gaussian integers and \ak — ak\ < 2“^ for all k = 0, ...,n. We 
also say that dk approximates to L bits after the binary point, and a 
corresponding polynomial F = do + ■ ■ ■ + dn ■ with coefficients fulfilling 
the latter properties is called an (absolute) L-bit approximation of F. It is 
assumed that the cost for asking the oracle for such an approximation is 
the cost for reading the approximations. 

• For any non-negative integer k, we denote by [A:] the set {1... A:} of size k. 
For any set S and any non-negative integer k, we write (^) for the set of 
all subsets of S of size k. 

• maxi(a;i,..., a;/j) := max(l, jcci |,..., |xfe|) for arbitrary Xi,...,Xk G C, 
log := log 2 the binary logarithm, and 


log(a:i,... ,Xfc) := |■maxl(logmaxl(a:l,.. ■,Xk))^■ 


Notice that, if \z\ < 2 for some z G C, then log(z) is 1. Otherwise, log( 2 ;) 
equals log \z\ rounded up to the next integer. 

• Halloo = max{|afe| : k = 0, ...,n} denotes the infinity-norm of F. We 
further define Tp := ldg(||F||oo), which bounds the number of bits before 
the binary point in the binary representation of any coefficient of F. 
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• r^;’ := lc)g(max"_]^ \zi\) is defined as the logarithmic root bound of F. 

• Mea^ := \an\ ■ riti inaxi(2i) is defined as the Mahler measure of F. 

• ap{zi) := minj^i \zi — Zj\ is defined as the separation of the root Zi and 
ap ■= min"^]^ ap{zi) as the separation of F. 

• For an arbitrary region 7^ C C in the complex space, we define <jp(TZ) := 
mini-zieTZ^Fizi), which we call the separation of F restricted to TZ. We 
further denote by Z{TZ) the set of all roots of F that are contained in TZ, 
and by Meai?(7^) := |a„| • Hz g 2(7 ?.)Mahler measure of F 
restricted to TZ. 

• We denote the interior of a disk in the complex plane with center m € C 
and radius r G M"*" by A = A(m,r). For short, we also write A • A to 
denote the disk A(m, A • r) that is centered at m and scaled by a factor 
A S K"*". We further use Fa (x) to denote the shifted and scaled polynomial 
F{m + r ■ x), that is, F^ix) := F{m + r ■ x). 

• A disk A is isolating for a root Zi of F if it contains Zi but no other root 
of F. For a set S of roots of F and positive real values pi and with 
Pi < 1 < /92, we further say that a disk A is (pi, P 2 )-isolating for S' if pi • A 
contains exactly the roots contained in S and p2 ■ A \ pi • A contains no 
root of F. 

• Throughout the paper, we only consider squares 

B = ^Z = X i • y G ^ X G [Xniin, ^max] and y G [pmin, ymax]} 

in the complex space that are closed, axis-aligned, and of width w{B) = 2^ 
for some £gZ (i.e., jx^ax - a^minl = l2/max - 2/min| = 2^), hence, for brevity, 
these properties are not peculiarly mentioned. Similar as for disks, for a 
A G M"*", A • B denotes the scaled square of size A • 2^ centered at B. 

According to Cauchy’s root bound (e.g. see [62]), we have \zi\ < 1 + 
max'Fg < 1 + 4 • 2’’’^, and thus F 1 ? = 0{tp). In addition, it holds that 

Tp < lbg(2" • Meai?) < n(l + Fi?) < 2nr p. 

Following [28, Theorem 1] (or [44, Section 6.1]), we can compute an integer 
approximation fi? G N of Fp with T^ + l < fi? < Fj? + 8 logn +1 using 0{n^Tp) 
many bit operations. For this, the coefficients of F need to be approximated to 
©(nrp) bits after the binary point. From F p, we then immediately derive an 
integer F = 2''', with 7 := [logf j?] G N>i, such that 

ri7' + i<f'p’<r < 2 'f'i?< 2'(ri? + 8 log u +1). (2) 

It follows that 2^ = is an upper bound for the modulus of all roots 

of F, and thus once can always restrict the search for roots to the set of all 
complex numbers of absolute value of at most 2^. 
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3 Counting Roots in a Disk 

In this section, we introduce the T*(A)-test, which constitute our main ingredient 
to count the numbers of roots of in a given disk A. Here, we briefly summarize 
the main properties of the T*(A)-test. The reader willing to focus on the 
algorithmic details of the root isolation algorithm is invited to read the following 
summary and skip the remainder of this section on a first read. 

• For a given polynomial F as in (1) and a disk A, the T*(A)-test always 
returns an integer k G { — 1,0,1,, n}. If fc > 0, then A contains exactly 
k roots of F. If A: = — 1, no further information on the number of roots in 
A can be derived; see Lemma 4, part (b). 

• If A is (pi,/92)-isolating for a set of k roots of F, where pi = « 0.94 

and P2 = |) then T*(A) returns k, see Lemma 4, part (a). In particular, 
T*(A) returns 0 if | • A contains no root. 


• The cost for the T* (A)-test is bounded by 

d(n(Ti. + nlog(m,r) +lbg(||FA||;;o^))) 

= d{n(TF +nlog(m,r) + lbg((max |F(z)|)“^))) 


zGA 


bit operations, and thus directly related to the size of A and the maximum 
absolute value that F takes on the disk A. For this, the test requires an 
L-bit approximation of F, with 


L = 0(rF + nlog(m,r) +log(|jFA||oo^)) 

= 0(rF + nlog(m,r) + l6g((max |F(z)|)"^)), 

z^A 


see Lemma 5. Here, we used that max^gA |^"(^)| A + 1) ■ ||-Pa||oo as 
shown in (27) in the proof of Theorem 3. 

3.1 Pellet’s Theorem and the T^-Test 

In what follows, let k be an integer with 0 < k < n = deg F, and let AT S K with 
K > 1. We consider the following test, which allows us to compute the size of a 
cluster of roots contained in a disk A(m,r): 

Definition 1 (The Tfc-Test). For a polynomial F € C[a;], the T^-test on a disk 
A A{m,r) with parameter K holds if 



(3) 


or, equivalently, if F^^'^ (m) 0 and 



12 














Mostly, we will write Tk{A, K, F) for Tk{m,r,K,F), or simply Tk{A,K) if 
it is clear from the context which polynomial F is considered. Notice that if 
the Tfe-test succeeds for some parameter K = Kq, then it also succeeds for any 
K with K < Kq. Clearly, Tk{m,r, K, F) is equivalent to Tk{0,l, K, F/\), with 
Fa (a:) := F(m + r ■ x). 

The following result is a direct consequence of Pellet’s theorem, and, in our 
algorithm, it will turn out to be crucial in order to compute the size of a cluster 
of roots of F; see [39, Section 9.2] or [61] for a proof. 

Theorem 1. If Tk{m, r, K, F) holds for some K gM. with K > 1 and some k G 
{0,..., n}, then A{m, r) contains exactly k roots of F counted with multiplicities. 

We derive criteria on the locations of the roots zi,..., of F under which 

the Ffc-test is guaranteed to succeed: 

Theorem 2. Let k be an integer with 0 < k < n = deg(F), let K G ^ with 
K > 1, and let ci and C 2 be arbitrary real values fulfilling 


C 2 • n Tn 


1 + 2F 
2K 


> Cl ■ n> 


maxi(/c) 

ln(l+8 f)’ 


(5) 


For a disk A = A(m,r), suppose that there exists a real X with 


X > max(4c2 • maxi(fc) • n^, 16K ■ maxi(fc)^ • n) 

such that A is (1, X)-isolating for the roots zi,... ,Zk of F, then Tk{cin ■ A, K, F) 
holds. 

In our algorithm, we will only make use of Corollary 1, which is actually a 
consequence of Theorem 2 with the specific values F := |, Ci := 16, C 2 := 64, 
A = 256n^, and thus ~ 12.49 • maxi(fc) and In ~ 0.29. 

Corollary 1. Let A be a disk in the complex space that is 16n'^)-isolating 
for a set of k roots (counted with multiplicity) of F. Then, Tk{A, |,F) holds. 

The proof of Theorem 2 is given in the appendix. In the proof, we separately 
bound the two sums in (4). We also derive a bound on the minimal distance 
between a root of the fc-th derivative F^^^ of F and a cluster of k roots of F. 
Pawlowski [37] provides a similar but more general bound, which implies a bound 
on the first sum in (4). However, compared to [37], our proof is significantly 
shorter and uses only simple arguments, hence we decided to integrate it in the 
appendix of this paper for the sake of a self-contained presentation. 


3.2 The T^-Test: Using Graeffe Iteration 

Corollary 1 guarantees success of the Ffc(A, 3/2, F)-test, with k = |Z(A)|, if 
the disk A is 16n‘^)-isolating for a set of k roots. In this section, we use 
a well-known approach for squaring the roots of a polynomial, called Graeffe 
iteration [4], in order to improve upon the T^-test. More specihcally, we derive a 
variant of the F^-test, which we denote T^-test^, that allows us to exactly count 
the roots contained in some disk A if A is (pi,/ 02 )-isolating for a set of k roots, 
with constants pi and p 2 of size pi « 0.947 and p 2 = |- 

^The superscript “G” indicates the use of Graeffe iteration. 
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Algorithm 1: Graeffe Iteration 

Input : Polynomial F{x) = a non-negative integer N. 

Output : Polynomial (x) = If F has roots zi,... ,Zn, 

then has roots zf",, 0 ^”, and 

1 F[01(x) := F{x) 

2 for i = 1,..., A do 

3 1^ FW(a;) := (—— x ■ Fo~^\xY] 

4 return F^^\x) 


Definition 2 (Graeffe Iteration). For a polynomial F{x) = ^ *C[a;], 

write F{x) = Ff.{x^) -I- x ■ Fo{x‘^), with 

Fe(x) := a2i!i\x^^^ + -f ... -f a2X + oq, and 

Fo{x) := -|- -f ... -|- a^x + ai. 

Then, the first Graeffe iterate of F is defined as: 

FW(cr) := i-inF^ixf - X ■ F^ixn 

The first part of the following theorem is well-known (e.g. see [4]), and we 
give its proof only for the sake of a self-contained presentation. For the second 
part, we have not been able to find a corresponding result in the literature. 
Despite the fact that we consider the result to be of independent interest, we 
will need it in the analysis of our approach. 

Theorem 3. Denote the roots of F by zi, ..., Zn, then it holds that (x) = 

Er=oar'2;* = al-Utii X — zf). In particular, the roots of the first Graeffe 
iterate are the squares of the roots of F. In addition, we have 

n2.maxi(||F|U)2>||F^W||^>||F||L-2-4T 

Proof. See Appendix 7.2. □ 

We can now iteratively apply Graeffe iterations in order to square the roots of 
a polynomial F(x) several times. In this way, we can now reduce the “separation 
factor of the Tfe-Test” from polynomial in n (namely, 256n®) to a constant value 
(in our case, this constant will be « 1.41) when we run N, with N = ©(loglogn), 
Graeffe iterations first, and then apply the T^-test; see Algorithm 2. From 
Theorem 2 and Theorem 3, we then obtain the following result: 

Lemma 1. Let A be a disk in the complex plane and F{x) € C[a;] a polynomial 
of degree n. Let 

N |'log(l -I- logn)] -I- 5 ( 6 ) 

and 

Pi ■= « 0.943 and p 2 := ^ (7) 

Then, we have > pi, and it holds: 
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Algorithm 2: T^(A, A)-Test 

Input : Polynomial F{x) of degree n, disk A = A(m,r), real value K 
with 1 < A < I 

Output : True or False. If the algorithm returns True, A contains exactly 
k roots of F. 

1 Call Algorithm 1 with input Fa{x) := F{m + r ■ x) and 
N := [log(l + logn)] + 5, which returns F^^ 

2 return Tfe(0,1, A, (x)) 


(a) If A is (pi, P 2 )-isolating for a set of k roots of F, then T®(A, |) succeeds. 

(b) //T^(A, A) succeeds for some A > 1, then A contains exactly k roots. 


Proof. The lower bound on p{n) := ^ y jA follows by a straight forward compu¬ 
tation that shows that p{n), considered as a function in n, is strictly increasing 
and that p{2) « 0.947 > « 0.943. Now, let F^^ be the polynomial obtained 

from Aa after performing N recursive Graeffe iterations. If A is (pi, p 2 )-isolating 
for a set of k roots of F, then the unit disk A' := A(0,1) is also (pi, p 2 )-isolating 
for a set of k roots of Fa, that is. A' contains k roots of Fa and all other roots 
of Fa have absolute value larger than |. Hence, we conclude that F^^ has k 
roots of absolute value less than p^” < jF, whereas the remaining roots have 
absolute value larger than p^ > 16n"^. From Corollary 1, we thus conclude that 
Tfe(A', a]^^) succeeds. This shows (a). Part (b) is an immediate consequence 
of Theorem 1 and the fact that Graeffe iteration does not change the number of 
roots contained in the unit disk. □ 

In the special case where fc = 0, the failure of (A) already implies that 
I • A contains at least one root. 


The following result is a direct consequence of Theorem 3. We will later use 
it in the analysis of our algorithm: 

Corollary 2. Let Fa and A_^^ be defined as in Algorithm 2. Then, it holds: 
1^(I1-P'a^'(2;)||oo, II (a;) 11”^) = 0(logn • (n -b lbg(|| Aa||oo, ||-f^A||;;o^))- 


3.3 The T^-Test: Using Approximate Arithmetic 

So far, the T^-test is formulated in a way such that, in general, high-precision 
arithmetic, or even exact arithmetic, is needed in order to compute its output. 
Namely, if the two expressions on both sides of (3) are actually equal, then 
exact arithmetic is needed to decide equality. Notice that, in general, we cannot 
even handle this case as we have only access to (arbitrary good) approximations 
of the coefficients of the input polynomial A. But even if the two expression 
are different but almost equal, then we need to evaluate the polynomial A and 
its higher order derivatives with a very high precision in order to decide the 
inequality, which induces high computational costs. This is a typical problem 
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that appears in many algorithms, where a sign predicate V is used to draw 
conclusions, which in turn decide a branch of the algorithm. Suppose that, 
similar as for the T^-test (with E( = 

there exist two non-negative expressions Ei and E^ such that V succeeds® if and 
only if Ei — E^ has a positive sign (or, equivalently, \i Ei > E^). We further 
denote by the predicate that succeeds if and only if the stronger inequality 
Ei — ^ ■ Er >0 holds.® Then, success of Va implies success of V', however, 
a failure of Va does, in general, not imply that V fails as well. As already 
mentioned above for the special case, where V = Tk{m^r^l, E), it might be 
computationally expensive (or even infeasible) to determine the outcome of V, 
namely in the case where the two expressions Ei and E^ are equal or almost 
equal. In order to avoid such undesirable situations, we propose to replace the 
predicate 7^ by a corresponding so-called soft-predicate [61], which we denote 
by V. V does not only return True or False, but may also return a flag called 
“Undecided”. If it returns True or False, the result of V coincides with that 
of V. However, if V returns Undecided, we may only conclude that Ei is a 
relative |-approximation of (i.e., ^ ■ Ei < E^ < ^ ■ Ei). We briefly sketch 
our approach and give details in Algorithm 3: In the first step, we compute 
approximations Ei and E^ of the values Ei and E ^, respectively. Then, we check 
whether we can already compare the exact values Ei and Er by just considering 
their approximations and taking into account the quality of approximation. If 
this is the case, we are done as we can already determine the outcome of V. 
Hence, we define that V returns True (False) if we can show that Ei > E^ 
{Ei < Er). Otherwise, we iteratively increase the quality of approximation until 
we can either show that Ei > Er, Ei < Er, or ^ ■ Ei < Er < ^ ■ Ei. We may 
consider the latter case as an indicator that comparing Ei and Er is difficult, 
and thus V returns Undecided in this case. 

It is easy to see that Algorithm 3 terminates if and only if at least one of 
the two expressions Ei and Er is non-zero, hence we make this a requirement. 
In the following lemma, we further give a bound on the precision to which the 
expressions Ei and Er have to be approximated in order to guarantee termination 
of the algorithm. 

Lemma 2. Algorithm 3 terminates for an L that is upper bounded by 
Lo := 2 • {\^{ma.x{Ei, Er)~^) 4). 

Proof. Suppose that L > \og{max{Ei, Er)~^) 4. We further assume that 

Ei = rriax.{Ei, Er)] the case Er = max{Ei, Er) is then treated in analogous 
manner. It follows that 

E+ <Er + 2-^+1 <Ei + 2-^+1 <l.Ei<\-Ei- 2 "^+® <\-E^. 

o Z Z 

Hence, if, in addition, | • E'^ < Eif, then the algorithm returns Undecided in 
Step 10. Otherwise, we have ^ • Ei > Ei > E'^ > | • Ef, and thus 

Ef >Ei- 2-^+1 >-■ Ei>- ■ EiE 2-^+^ > Ef + 2"^+^ > E+, 

8 4 

^We assume that the predicate V either returns “True” or “False”. We say that V succeeds 
if it returns True. Otherwise, we say that it fails. 

®You may replace | by an arbitrary real constant K larger than 1. 
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Algorithm 3: Soft-predicate V 

Input : A predicate V defined by non-negative expressions Ei and E^, 
with Eg ^ 0 or Ej. ^ 0; i.e. V succeeds if and only if Eg > E^. 
Output : True, False, or Undecided. In case of True (False), V succeeds 
(fails). In case of Undecided, we have ■ Eg < Er < ^ ■ Eg. 


1 

2 

3 

4 

5 

6 

7 

8 


L ■= 1 

while True do 

Compute L-bit approximations Eg and Er of the expressions Eg and 
Er, respectively. 

Ef^ := max(0. Eg ± 2“^) and E^ := max(0, Er ± 2“^) 


if Ef > i?+ then 
1^ return True 

if E^ < Et then 
1^ return False 


// It follows that Eg > Er- 


9 

10 

11 


// It follows that Eg <Er. 
if I ■ E^ < Ef < E+ < ^ ■ Ef, then 
1^ return Undecided 

// It follows that ■ Eg < Er < ^ ■ Eg. 

L:=2- L 


which shows that the algorithm returns True in Step 6. Since we double L 
in each iteration, it follows that the algorithm must terminate for an L with 
L < 2 ■ {\og{ma.x{Eg, Er)~^) + 4). □ 

Notice that if V returns True, then V also succeeds. This however does 
not hold in the opposite direction. In addition, if Va succeeds, then Eg > Er 

and Eg cannot be a relative |-approximation of Er, hence V must return True. 
We conclude that our soft-predicate is somehow located “in between” the two 
predicates V and U|. 

We now return to the special case, where V = Tk{m,r, 1,F), with Eg = 

and Er = expressions on the left and 

the right side of (3), respectively. Then, success of V implies that the disk 
A = A(m,r) contains exactly k roots of F, whereas a failure of V yields no 
further information. Now, let us consider the corresponding soft predicate 
V — Tk{A,F) oi V = Tk{A,F). If V returns True, then this implies success 
of V. In addition, notice that success of Tk{A, |,T) implies that V returns 
True, and thus we may replace Tfc(A, 1,^) by T);(A,F) in the second part of 
Theorem 2. Similarly, in Lemma 1, we may also replace Tjf{A, |,T) by the 
soft-version T^{A,F) of T®(A,F). We give more details for the computation of 
Tk{A, F) and T® (A, F) in Algorithms 4 and 5, which are essentially applications 
of Algorithm 3 to the predicates Tk{A,F) and T^{A,F). The lemma below 
summarizes our results. Based on Lemma 2, we also provide a bound on the 
precision L for which Algorithm 4 terminates and a bound for the bit complexity 
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Algorithm 4: Tfc(A,F)-test 

Input : A polynomial F{x) of degree n, a disk A := A(to, r) in the 
complex plane, and an integer k with 0 < k < n. 

Output : True, False, Undecided. If the algorithm returns True, the disk 
A(m, r) contains exactly k roots of F. 


1 L:=l 

2 while True do 

Compute an approximation Fa (a;) = polynomial 

Fa(x) := F{m + r ■ x) such that fi ■ g Z 

and \fi- fi\< for all i. 

// [L + [log(n + l)~\)-bit approximation of F/^. 
f~ := max(0, \fi\ — ) for * = 0,..., n. 

/+ := l/^l + 2 --E^-ri°g("+i)l for i = 0,... ,n. 

// lower and upper bounds for \ fi\. 

if fk - E^^k /^ > 0 then 
1^ return True 

// It follows thatTk{A,F) succeeds. 

if E.^kfr-fk 

1^ return False 

// It follows that Tj;(A,U) fails. 

10 if T,^^k fr - I ■ fk (^rid I • /- - /^ > 0 then 

11 1^ return False 

12 L := 2 ■ L 


fjf>0 then 


of Algorithm 4. A corresponding bound for the bit complexity of carrying out 
the T^{A, F)-test for all /c = 0,... , n is given in Lemma 3. 

Lemma 3. For a disk A := A{m,r) in the complex plane and a polynomial 
F g C[a;] of degree n, the Tk{A, F)-test terminates with an absolute precision L 
that is upper bounded by 

L{A, F) := L{m, r, U) := 2 ■ (4 + lbg(||FA||-i)) . (8) 

If Tk{A, F) succeeds, the Tk{A, F)-test returns True. The cost for running 
the Tk{A, F)-test for all fc = 0,..., n is upper bounded by 

0(n(n • lbg(m, r) + tf + L{A, F))) 

bit operations. The algorithm needs an 0{n ■ log{m,r) + + L{A, F))-bit 

approximation of F. 

Proof. Let V := Tk{A,l, F) he the predicate that succeeds if and only if Eg > E^, 
with El := \fk\ and Er := Y.i^k \M- Then, Ef := and Ef^ := f^ are 

lower and upper bounds for Eg and Er, respectively, such that \E^ — Ei\ < 2“^+i 
and \Ef — Er\ < 2“^+i. Hence, Lemma 2 yields that Algorithm 4 terminates 
for an L smaller than 2 • (4 + lbg(max(Fif, Er)~^)) < L{A, F). 
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Algorithm 5: T^(A, F)-Test 

Input : Polynomial F{x) G C[a;] of degree n, a disk A := A(m, r) in the 
complex space. 

Output : True, False, or Undecided. If the algorithm returns True, A 
contains exactly k roots of F. 

1 Let be the N-th Graeffe iterate of Fa (a;) := F{m + r • cc), where 

N := |'log(l + logn)] + 5 

2 Output ffe (0,1, 


We have already argued above that success of the predicate Fa = Ffe(A, |, F) 
implies that P = Ffc(A, F) returns True. Hence, it remains to show the claim on 
the bit complexity for carrying out the Tk{A, F)-test for all fc = 0,..., n. For 
a given L, we can compute an (F + [log(n + l)])-bit approximation Fa (a;) = 
Sr=o with a number of bit operations that is bounded by 0{n{TF + 

n lbg(TO, r) + L)); e.g. see the first part of the proof of [45, Lemma 17]. For a fixed 
k, the computation of the signs of the sums in each of the three IF clauses needs 
n additions of dyadic numbers with denominators of bit size [log(n + 1)] + L 
and with numerators of bit size 0{L + nlbg(r) + Tp), hence the cost is bounded 
by 0{n(Tp + nlbg(r) + L)) bit operations. Notice that, when passing from an 
integer k to & k' ^ k, the corresponding sums in one IF clause differ only by two 
terms, that is, and /^. Hence, we can decide all IF clauses for all k using 
0(n) additions. Furthermore, we double the precision L in each step, and the 
algorithm terminates for an L smaller than F(A,F). Hence, L is doubled at 
most logF(A, F) many times, and thus the total cost for all k is bounded by 
0{n{Tp + n lbg(TO, r) + F(A, F))) bit operations. □ 

We now extend the above soft-variant of the T^-test to a corresponding 
soft-variant of the T^-test, which we denote T®; see Algorithm 5 for details. We 
further combine for all A: = 0,..., n to obtain T*(A, F) with 

\ k if there exists a k such that T^{A) succeeds 

T*(A, F) ■= \ , . (9) 

I —1 otherwise. 

Again, for brevity, we often omit F and just write T*(A). We say that T* 
succeeds if it returns a non-negative value. Otherwise, it fails. 

The following result, which can be considered as the “soft variant” of Lemma 1, 
can then immediately be deduced from Lemma 1 and Lemma 3: 

Lemma 4 (Soft-version of Lemma 1). Let A := A(m, r) be a disk in the complex 
plane, F{x) G C[a;] be a polynomial of degree n, and let pi = and P 2 = |- 
Then, it holds: 

(a) If A is {pi, P 2 )-isolating for a set of k roots of F, then T»(A) returns k. 

(b) //T*(A) returns a k>0, then A contains exactly k roots. 

For the complexity analysis of our root isolation algorithm (see Section 4) , 
we provide a bound on the total cost for running the T*-test. 


19 





Lemma 5. The total cost for carrying out the T*(A) is bounded by 


d{n{Tp + n\og{m,r) + L{A, F))) = 0{n{Tp +nlog(m, r) + log((max |F(z)|) ^))) 

A 

bit operations. For this, we need an L-bit approximation of F with 
L = d{Tp + n\og{m,r) + L{A, F)) = 0(tf+ nlog(m, r)+log((max |i^( 2 :)|)“^)). 

zGA 

Proof. According to Lemma 3, the computation of Tfe(0,1, F^^) needs an L-bit 
approximation F^^ of F^^ , with L bounded by 

d{n + + L(0,1, Lf')) = 0(n + lbg(||Ff > ||oo, H^f' ||-')). (10) 

Given such an approximation F^^, the cost for running the test for all k = 

0,..., n is then bounded by 0(n{n + rp[iv] -I- L)) bit operations. In each of the 

N = O(loglogn) Graeffe iterations, the size of IbgdlFj^^Hoo, ||.Fa lloo^) increases 
by at most a factor of two plus an additive term 4n; see Theorem 3. Hence, we 
must have 

logdl-f^A ||oo,||^a'||;)o^) = 0(logn •ldgdiLA||oo,||iA||;;o^) +nlogn) 

= 0(n ldg(TO, r) + Tp + L(A, F)) 

for all i = 0,...,N. We conclude that the above bound (10) for L can be 
replaced by 0{tp F nlbg(m,r) -|- L{A,F)). 

It remains to bound the cost for computing an approximation of f]^^ 
with II— Fa Iloo < 2“^. Suppose that, for a given p G N we have computed 
an approximation Fa of Fa, with ||Fa — Fa||oo < 2 “'’. According to [47, 
Theorem 8.4] (see also [23, Theorem 14] and ]45, Lemma 17]), this can be achieved 
using a number of bit operations bounded by 0(n(n lbg(m, r) + Tp + p)). In 
each Graeffe iteration, an approximation F^^ of f](^ is split into two polynomials 
F^o ^A e coefficients of comparable bit size (and half the degree), 
and an approximation of f](^ is then computed as the difference of f'^^ 

and X ■ F ^If all computations are carried out with fixed point arithmetic 
and an absolute precision of p bits after the binary point, then the precision 
loss in the Fth step, with t = 0,..., A, is bounded by 0(log n + log \\F^ ||oo) = 
0(2®(logn + log ||Fa||oo)) = 0(logn(logn + log ||Fa||oo)) bits after the binary 
point. The cost for the two multiplications and the addition is bounded by 
0{n{p + log ||f](^||oo)). Since there are only N = O(loglogn) many iterations, 
we conclude that it suffices to start with an approximation Fa of Fa, with 
||Fa — FaIIoo < 2“^ and p = 0(nlog(m,r) + -|- L(A,F)). The total cost for 

all Graeffe iterations is then bounded by 0{np) bit operations, hence the claim 
follows together with the fact that max^gA l^(■z)| ^ {n + 1)||Fa||oo as shown 
in (27) in the proof of Theorem 3. □ 
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4 CISOLATE: An Algorithm for Root Isolation 

We can now formulate our algorithm, which we denote by CIsolate, to isolate 
all complex roots of a polynomial F{x) that are contained in some given square^® 
B C C. If the enlarged square 2B contains only simple roots of F, then our 
algorithm returns isolating disks for all roots that are contained in B. However, 
it might also return isolating disks for some of the roots that are not contained 
in B but in the complement 2B \ B. In particular, in the important special case, 
where F is square-free and where we start with a square B that is known to 
contain all complex roots of F, our algorithm isolates all complex roots of F. 
Before we give details, we need some further definitions, which we provide in 
Section 4.1. In Section 4.2, we first give an overview of our algorithm before we 
provide details and the proof for termination and correctness. 

4.1 Connected Components 

Given a set S = {Bi ,..., Bm} of squares Bi,..., Bm C C, we say that two 
squares B,B' G S are connected in S {B ^5 B' for short) if there exist squares 
Bi-^,..., Bi^, G S with Bi-^ = B, Bi^, = B', and Bi. n Bi.^^ ^ 0 for all 
j = 1,... ,s' — 1. This yields a decomposition of S into equivalence classes 
Cl,... ,Ck G S that correspond to maximal connected and disjoint components 
Ce = Ui-s eCf £ = 1,... ,k. Notice that Cg is defined as the set of 

squares Bi that belong to the same equivalence class, whereas Ce denotes the 
closed region in C that consists of all points that are contained in a square 
Bi G Ci- However, for simplicity, we abuse notation and simply use C to denote 
the set of squares B contained in a component C as well as to denote the set 
of points contained in the closed region C. Now, let C = {Bi,... ,Bs] be a 
connected component consisting of equally sized squares Bi of width w, then we 
define (see also Figure 1): 

• Be is the axis-aligned closed square in C of minimal width such that 
C C Be and 

min 5ft(z) = minSftfz) and max 9(z) = max3(z), 
z£Bc zdC zdBc zdC 

where 5R(z) denotes the real part and 3 ( 2 :) the imaginary part of an 
arbitrary complex value z. We further denote me the center of Be, and 
Ac := A(mc, \w{Be)) a disk containing Be, and thus also C. 

We further define the diameter w{C) of the component C to be the width 
of Be, i.e. w{C) := w{Be), and r{C) := to be the radius of C. 

• C~^ := 2i3i is defined as the union of the enlarged squares 2Bi. 

Notice that is the ^-neighborhood of C (w.r.t. max-norm). 

4.2 The Algorithm 

We start with an informal description of our algorithm CIsolate, where we 
focus on the main ideas explaining the ratio behind our choices. For the sake 

'^'^As already mentioned in Section 2, we only consider closed, axis-aligned squares B G C. 
Hence, these properties are not further mentioned throughout the following considerations. 
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of comprehensibility, we slightly simplified some steps at the cost of complete 
formal correctness, hence, the considerations below should be taken with a grain 
of salt. A precise definition of the algorithm including all details is given in 
Algorithm 6 and the subroutines NewtonTest (Algorithm 7) and Bisection 
(Algorithm 8). 

From a high-level perspective, our algorithm follows the classical subdivision 
approach of Weyl [57]. That is, starting from the input square B, we recursively 
subdivide B into smaller squares, and we remove squares for which we can show 
that they do not contain a root of F. Eventually, the algorithm returns regions 
that are isolating for a root of F. In order to discard a square B, with B C B, 
we call the T*(Ab,F)- test^^, with the disk containing B. The remaining 
squares are then clustered into maximal connected components. We further check 
whether a component C is well separated from all other components, that is, we 
test whether the distance from C to all other components is considerably larger 
than its diameter. If this is the case, we use the T*-test in order to determine 
the “multiplicity” kc of the component C, that is, the number of roots contained 
in the enclosing disk Ac; see Line 9 of Algorithm 6 and Figure I for details. 
If fee = 1, we may return an isolating disk for the corresponding unique root. 
Otherwise, there is a cluster consisting of two or more roots, which still have to 
be separated from each other. A straight-forward approach to separate these 
roots from each other is to recursively subdivide each square into four equally 
sized squares and to remove squares until, eventually, each of the remaining 
components contains exactly one root that is well separated from all other roots; 
see also Algorithm 8 (Bisection) and Figure 3. However, this approach itself 
yields only linear convergence to the roots, and, as a consequence, there might 
exist (e.g. for Mignotte polynomials) long sequences Ci,... ,Cs of interlaced 
connected components with invariant multiplicity A:, that is Ci D (^2 D • ■ ■ D 
and k = kc^ = • • • = fee, > 1- The main idea to traverse such sequences more 
efficiently is to consider a cluster of k roots as a single root of multiplicity k and 
to use Newton iteration (for multiple roots) to compute a better approximation 
of this root. For this, we use an adaptive trial and error approach similar to the 
quadratic interval refinement (QIR) method, first introduced by Abbott [1]; see 
Algorithm 7 (NewtonTest) and Figure 2. In its original form, QIR has been 
combined with the secant method to efficiently refine an interval that is already 
known to be isolating for a real root of a real polynomial. Recent work [42] 
considers a modified approach of the QIR method that uses Newton iteration (for 
multiple roots) and Descartes’ Rule of Signs. It has been refined and integrated 
in almost optimal methods [43, 45] for isolating and approximating the real 
roots of a real (sparse) polynomial, where it constitutes the crucial ingredient 
for quadratic convergence. In this paper, we further extend the QIR approach 
for approximating complex roots of a polynomial. 

The main crux of the NewtonTest (and the QIR method in general) is that 
we never have to check in advance whether Newton iteration actually yields an 
improved approximation of the cluster of roots. Instead, correctness is verified 
independently using the T»-test. In order to achieve quadratic convergence 
in the presence of a well isolated root cluster, we assign, in each iteration, an 
integer Nc to each component C. The reader may think of Nc as the actual 

fact, it suffices to just call the {Ab , F)-test. However, since and T+ have 
comparable complexity, we just stick to T* to simplify the presentation. 
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Algorithm 6: CIsolate 

Input : A polynomial F{x) S C[a;] as in (1) and a square S C C of width 
Wq := w{B) = 2^°, with £q G Z-, F has only simple roots in 2B. 

Output : A list O of disjoint disks Ai,..., Ag C C such that, for each 
f = 1,..., s, the disk A^ as well as the enlarged disk 2Ai is 
isolating for a root of F that is contained in 2B. In addition, for 
each root z G B, there exists a disk Ai G O that isolates z. 


iO = {} 

2 C = {iBA)} 


// list of isolating disks 

// list of pairs {C,Nc), with C a connected com- 
/ / ponent eonsisting of sc equally sized squares, 
/ / each of width 2^^', where £c € ^<4 ■ 

// integer with Nq = 2^ and nc G N>i. 


// * Preprocessing * / / 

3 repeat 


Let {C, Nc) be the unique pair in C 

// V [Jc-{c Nc}&C C = B, then there exists 
// unique component C with {C,Nc) G C. 


// * linear step *// 

5 {C{,...,C'A := Bisection(C') and C = {(Ci, 4),..., (C;, 4)} 

6 until Uc:(c,Arc)eC ^ ^ B 


// * Main Loop *// 

7 while C is non-empty do 


9 


10 

11 

12 


13 


14 


Remove a pair [C, Nc) from C. 

if 4Ac n C" = 0 for eaeh (C", Nc') G C with C ^ C and there exists a 
kc G {I,... ,n} such that kc = T*(2Ac’) = T*(4Ac') then 

// If the second condition holds, kc equals the 
// number of roots eontained in 2Ac and 4Ac<. 

if fcc = 1 then 

1^ Add the disk 2Ac to O, continue 
if Ac > 1 then 

Let xc G B\C he an arbitrary point with distance 2^*^“^ 
from C and distance 2^'^“^ or more from the boundary of B. 

// Existence of such a point follows from the 
// proof of Theorem f. It holds that F{xc) ^ 0. 
if NewtonTest(( 7, Ac, Ac,a:c) = (Success, C") then 

// * quadratic step *// 

Add {C',Nq) to C, continue 


// * linear step *// 

{C[,...,C'f}:= Bisection(C'). 

16 |_ Add (C^, max(4, \/Nc)), ..., {C[, max(4, \/Nc)) to C. 

17 return O. 


speed of convergence to the cluster of roots contained in C. Then, in case of 
success of the NewtonTest, the component C is replaced by a component 
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Algorithm 7: NewtonTest 

Input : A tuple (C, Nq, kc,xc)- C = {Bi,..., Bg^} is a connected 
component consisting of equally sized and aligned squares B^ 
contained in B and of size 2^'^, Nc is an integer of the form 2^ 
with nc G N>i, kc is the number of roots in 4Ac, and xc is a 
point with F{xc) ^ 0. 

Output : Either Failure or (SUCCESS, C"), where C" C C is a connected 
component that contains all roots contained in C. C consists of 
at most 4 equally sized and aligned squares, each of width ■ 


1 

2 

3 


if Algorithm 3 does not return False for the input Ef, := 4r(C)|-F''(a;c)| 
and Er := |F(a;c)| theu // This implies |_F(a;c)| < ^f'{C)\F'{xc)\- 

for L = 1, 2,4,... do 

Compute L-bit approximations of Fixe) and F'[xc) and derive an 
(6 — £c + log Ac)-bit approximation x'q of the Newton iterate 


x'c ■■= xc - kc- 


F{xc) 

F'{xc) 


1 2^c 

such that \^c - xc\ < - (11) 


4 

5 

6 

7 

8 

9 

10 

11 


// For more details, consider the similar 
// computation in [45, Step 2 of NewtonTest]. 

Let A':=A(Fc,l-g). 
if A' n C = 0 theu 
1^ returu Failure 

if T*(A') = kc holds 

// Then, A' contains all roots contained in 2Ac- 

theu 

Decompose each square Bi into 4A^ many equally sized 
sub-squares Bij 

returu (SUCCESS, C"), with C the unique connected component 
consisting of all squares Bij of width that intersect A'. 


12 returu Failure 


C' C C oi diameter w{C') « w{C) ■ . In this case, we “square the speed’ 

of convergence", that is, we set Nc ■= Nq. If the NewtonTest fails, we 
fall back to bisection and decrease the speed of convergence, that is, we set 
Nc' ■= '/Nc for all components C' into which the component C is split. Our 
analysis shows that the NewtonTest is the crucial ingredient for quadratic 
convergence. More precisely, we prove that, in the worst-case, the number s of 
components in each sequence Ci,..., as above becomes logarithmic in the 
length of such a sequence if only bisection would be used; see Lemma 8. 

We now turn to the proof of termination and correctness of the algorithm. 
In addition, we derive further properties, which will turn out to be useful in the 
analysis. 

Theorem 4. The algorithm CIsolate terminates and returns a correct result. 
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Algorithm 8: Bisection 

Input : A connected component C = {Bi, ..., Bg^} consisting of aligned 
squares Bi, each of width w{Bi) = 2^"^. 

Output: A list of components C' C C, each consisting of aligned and 
equally sized squares of width The union of all Cj 

contains all roots of F that are contained in C. 


1 

2 

3 


C ":=0 

for each Bi G C do 

Remove Bi from C and subdivide Bi into four equally sized 
sub-squares Bij, with j = 1,..., 4, and add these to C. 


4 

5 

6 
7 


for each B G C do 

ifT4AB)=0 // 

then 

1^ Remove B from C. 


This implies that B contains no root. 


8 Compute maximal connected components C[,.. .C'^ from the squares in C. 

9 return C[,.. .C[ 


In addition, at any stage of the algorithm, it holds that: 

(a) For any {C,Nc) G C, the connected component C consists of disjoint, 
aligned, and equally-sized squares Bi,... ,Bs^, each of width 2^"^ with 
some £c G 

(b) For any two distinct pairs {Ci,Nci) G C and ((72, G C, the distance 
between Ci and C 2 is at least max(2^‘^i, 2 ^c' 2 ). In particular, the enlarged 
regions and (7^ are disjoint. 

(c) The union of all connected components C covers all roots of F contained 
in B. In mathematical terms, 

F{z) ^0 for all z G B\ |J (7. 

C-.{C,Nc)GC 

(d) For each square B produced by the algorithm that is not equal to the initial 
square B, the enlarged square 2B contains at least one root of F. 

(e) Each component C considered by the algorithm consists of sc 

squares. The total number of squares in all components C is at most 
9-times the number of roots contained in 2B, that 

^ sc<9-\Zi2B)\. 

C:a(C,Arc)GC 

(f) Let {C,Nc) be a pair produced by the algorithm. The sequence of ancestors 
of a component C produced by the algorithm is recursively defined as 
follows. It consists of the component C from which C resulted followed 

will later prove that even the total number of squares produced by the algorithm in 
all iterations is near-linear in the number Z{2B) of roots contained in 2B. 
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Figure 1: A component Ci := C consisting of 5 squares i3i,... i? 5 , the enclosing 
square Be with center m := me and the disks Ac, 2Ac and 4Ac- The disk 
4Ac intersects the component C 2 but does not intersect the component C^. 

by the ancestors of C'. We denote with anc*(C'), the first ancestor of 
C for which the Newton Test was successful, if such exists, otherwise 
anc*(C') = B. Then w{C) < A Moreover, 

( 9 ) 2 ^xSmb) - 4 < Ac < crF(2B) := 

mini. 2 ;.g 2 B is the separation of F restricted to 2B. 


Proof. Part (a) follows almost immediately via induction. Namely, a component 
C consisting of squares of size 2^*^ is either replaced by a single connected 
component consisting of (at most 4) squares of width jNe in line 14 after 
NewtonTest was called, or it is replaced by a set of connected components 
C C C, each consisting of squares of size 2^^~^ in line 16 after Bisection was 
called. 

For (b), we can also use induction on the number of iterations. Suppose first 
that a component C is obtained from processing a component D in line 16. If 
C is the only connected component obtained from D, then, by the induction 
hypotheses, it follows that the distance to all other components C, with C"nll = 
0, is at least max(2^^, 2 ^c'') > rnax(2^'^, 2 ^cr'). If £) splits into several components 
Cl,..., Cs, with s > 1, their distance to any component C', with C' n Zl = 0, 
is at least max(2^^, 2^^') > max(2^‘^i, 2^^') for all i. In addition, the pairwise 
distance of two disjoint components Ci and Cj is at least 2^'^“^ = 2^'^i for 
all i. Finally, suppose that, in line 14, we replace a component H by a single 
component C. In this case, C C D and C consists of squares of width 2^~^/Ne- 
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Figure 2: The NewtonTest: If T*(A') = kc, with A' := A(ip, ^J' ), then 
A' contains exactly kc roots of F. Since T*(2Ac) = kc and (7+ C 2Ac, it 
follows that A' contains all roots contained in C'*'. The sub-squares Bij of width 
that intersect A' yield a connected component C of width at most 
2^^ /Nc < w{C)/Nc- In addition, all roots that are contained in C are also 
contained in C". Further notice that if x'q is contained in C, then A' intersects 
at most four squares Bij. Otherwise, it intersects at most three squares. In 
each case, the squares are connected with each other, and the corresponding 
connected component C" has width at most 2^’^ jNc < w{C)INc- 


Hence, the distance from C to any other component C is also lower bounded by 

max(2^c,2^c'). 

For (c), notice that in line 7 of Bisection, we discard a square B only if 
the T*(Ab) = 0. Hence, in this case, B contains no root of F. It remains to 
show that each root of F contained in C is also contained in C", where C' C C 
is a connected component as produced in line 14 after NewtonTest was called. 
If T*(A') = kc, then A' contains kc roots; see Lemma 4. Hence, since A' is 
contained in 2Ac, and since 2Ac also contains kc roots (as T*(2Ac) = kc 
holds), it follows that A' contains all roots that are contained in C. The disk 
A' intersects no other component C' ^ C as the distance from C to C is larger 
than 2^"^, and thus, by induction, we conclude that (A' n H) \ C contains no root 
of F. This shows that C already contains all roots contained in C. 

We can now prove part (d) and part (e). Any square B ^ B that is considered 
by the algorithm either results from the Bisection or from the NewtonTest 
routine. If a square B results from the Bisection routine, then the disk 
Ab = A{mB,w{B)) contains at least one root of F, and thus also 2B contains 
at least one root. If a square B results from the NewtonTest routine, then 
2B even contains two roots or more. Namely, in this case, T*(A') = kc holds 
for the disk A' = A(m',r'), with r' = ^w{B) and kc > 1, and thus A' contains 
kc roots. Since 2B contains the latter disk, 2B must contain at least kc roots. 
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Ac := A(m,^w(13c)) 



Figure 3: The Bisection routine: The green (brighter) sub-squares are all 
squares B for which T*(Ab) ^ 0. They are grouped together into three maximal 
connected components, which contain all roots contained in C. All other sub¬ 
squares are discarded. 


This shows (d). From (d), we immediately conclude that, for each component 
C ^ B produced by the algorithm, the enlarged component (7+ contains at least 
one root of F. In addition, since C~^ is contained in 2B, each of these roots must 
be contained in 2B. The first part in (e) now follows from the fact that, for a 
fixed root of F, there can be at most 9 different squares B of the same size such 
that 2B contains this root. From part (b), it follows that, for any two distinct 
components Ci and C 2 , the enlarged components and do not intersect, 
and thus the total number of squares in all components is upper bounded by 
9 • \Z{2B)\, which proves the second part in (e). 

For (f), we may assume that Nc > 4 as, otherwise, the inequality becomes 
trivial. Denote with C'i,...,Cs the sequence of ancestors of C, with Ci := 
anc*(C'), Cg = C, and Ci D Q+i. By definition of anc*(C), it holds that 
w(C 2 ) < since the step from Ci to C 2 is a quadratic step. It 

follows that 


w{C) = w{Cg) < < w{C2) < ^ ^ 

since Nc^ = for i = 2,..., s. 

We can now show that the algorithm terminates; the inequalities in (g) will 
then follow from the proof of termination: Suppose that the algorithm produces 
a sequence Ci, C 2 ,..., of connected components, with s > logn -|- 6 and 
Cl D (72 D • • • D Cs- If, for at least one index i G {1,..., s — 1}, Ci+i is obtained 
from Ci via a quadratic step, then w{Ci+i) < w{Ci)/Nci < w(Ci)/4. Hence, 
in this case, we also have w{Ca) < Now, suppose that each C^+i is 

obtained from Ci via a linear step, then each square in Ci has size and 

thus w{Cs) < 9n ■ 2^'^!“®+^ < . This shows that, after at most logn -I- 6 

iterations, the width of each connected component is halved. Hence, in order to 
prove termination of the algorithm, it suffices to prove that each component C 
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of small enough width is terminal, that is C is replaced by an isolating disk in 
line 10 or discarded in NewtonTest or Bisection. The following argument 
shows that each component C of width smaller than w that is 

not discarded is replaced by an isolating disk. We have already shown that (7+ 
must contain a root ^ of F, and thus we have \mc — ^| < 2w{C) < (Tj’(^)/ 16 
and rc < We conclude that the disks 2Ac and A{mc,8rc) are both 

isolating for Then, Lemma 4 guarantees that T*(2Ac) = 1 and T*(4Ac<) = 1 
hold. Hence, if 4Ac intersects no other component C" ^ C, then the algorithm 
replaces C by the isolating disk 2Ac in line 10, because the if-clause in line 9 
succeeds with kc = ^- It remains to show that the latter assumption is always 
fulfilled. Namely, suppose that AAq intersects a component C ^ C, and let 
B and B' be arbitrary squares contained in C and C", respectively. Then, the 
enlarged squares 2B and 2B' contain roots ^ and respectively, and ^ and 
must be distinct as C~^ and are disjoint. Hence, the distance between B 

and B', and thus also the distance S between C and C, must be larger than 
aF{2B) - 2^c - 2^0 = 2,2w - 2^^ - 2 ^c' > Slu, _ 2 ^c'. Hence, if 2 ^c' < 25ry, 
then 4Ac C A{mc,6w) does not intersect C'. Vice versa, if 2^c' > 25w, then 
the distance between C and C’ is at least max(2^‘^, 2^c') > 2bw, and thus 4Ac 
does not intersect C as well. Notice that (g) now follows almost directly from 
the above considerations. Indeed, \et C ^ B be an arbitrary component C 
and let D be any component that contains C. Since D is not terminal, we 
conclude that w{D) > w, and thus Nf> < according to (f). Since Nq 

is smaller than or equal to the square of the maximum of all values Nf, the 
second inequality in (g) follows. The first inequality follows from the fact that 

> minD-.CcD ^ 

For correctness, we remark that each disk D returned by the algorithm is 
actually isolating for a root of F contained in 2B and that 2D also isolates this 
root. Namely, for each component C produced by the algorithm, the enlarged 
component C+ contains at least one root. Now, if the if-clause in line 9 succeeds 
on C with kc = 1, it holds that T*(2Ac) = 1, and thus the disk 2Ac contains 
exactly one root Hence, since Ac contains C~^, this root must be contained 
in . In addition, if also T*(4Ac) = 1 holds, then the disk 4Ac< isolates ^ as 
well. Finally, it remains to show that the algorithm returns an isolating disk for 
each root ^ that is contained in B. From (a) and (c), we conclude that there is 
a unique maximal sequence S = Ci^C^, ■ ■ ■ ,Cs of connected components, with 
Cl D (72 D • • • D Cs, such that each Ci contains Now, when processing Cg, 
Cs cannot be replaced by other connected components C C Cg as one of these 
components would contain and this would contradict the assumption that the 
sequence S is maximal. Since Cg contains it cannot be discarded in Bisection 
or NewtonTest, hence Cg is replaced by an isolating disk for ^ in line 10. □ 

Remarks. We remark that our requirement on the input polynomial F to have 
only simple roots in 2B is only needed for the termination of the algorithm. 
Running the algorithm on an arbitrary polynomial (possibly having multiple 
roots) yields isolating disks for the simple roots as well as arbitrarily small 
connected components converging against the multiple roots of F in B. Namely, 
if B is not discarded in the first iteration, then the enlargement C'^ of each 
component C contains at least one root. Since C consists of at most 9n squares, 
each of size 2^'^ , it holds that each point in C approximates a root of F to an 
error of less than n ■ 2^'^+'^. In addition, the union of all components covers all 
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roots contained in B, and thus our algorithm yields L-bit approximations of 
all roots in B if we iterate until (-c < ~4 — logn — L for all components C. In 
the special situation, where we run the algorithm on an input square that is 
known to contain all roots and if, in addition, the number k of distinct roots of 
F is given as input, our algorithm can be used to return isolating regions for 
all roots. Namely, in this situation, we may proceed until the total number of 
connected components C equals k. Then, each of the enlarged components 
isolates a root of F. The latter problem is of special interest in the context of 
computing a cylindrical algebraic decomposition, where we have to isolate the 
roots of a not necessarily square-free polynomial with algebraic coefficients. In 
this case, it might be easier to first compute k via a symbolic pre-computation 
and to consider sufficiently good approximations of the initial polynomial instead 
of computing approximations of the square-free part of F. A corresponding 
approach based on approximate polynomial factorization has been presented 
in [28] , and we refer the reader to this work for more details and for a motivation 
of the problem. 

5 Complexity Analysis 

We split the analysis of our algorithm into two parts. In the first part, we focus 
on the number of iterations that are needed to isolate the roots of F{x) that are 
contained in a given square B. We will see that this number is near-linear^^ in 
\Z{2B)\, the number of roots contained in the enlarged square 2B. We further 
remark that, for any fixed non-negative constant e, the total number of iterations 
is near-linear in \Z{(1 -|- e) • B)\\ however, for the sake of simplifying analysis, 
we only provide details for the special case e = 1. Hence, we conclude that our 
algorithms performs near-optimal with respect to the number of subdivision 
steps if the input square B has the property that each root contained in -B 

is also contained in B] in particular, this is trivially fulfilled if B is chosen large 
enough to contain all roots of F. 

In the second part of our analysis, we give bounds on the number of bit 
operations that are needed to process a component C. This eventually yields a 
bound on the overall bit complexity that is stated in terms of the degree of F, the 
absolute values and the separations of the roots in Z{2B), and the absolute value 
of the derivative F' at these roots. For the special case, where our algorithm is 
used to isolate all roots of a polynomial of degree n with integer coefficients of 
bit size less than r, the bound on the bit complexity simplifies to 0{n^ -I- 

5.1 Size of the Subdivision Tree 

We consider the subdivision tree Tg, or simply T, induced by our algorithm, 
where B is the initial square/component. More specifically, the nodes of the 
(undirected) graph T are the pairs (C, Nc) G C produced by the algorithm, 
and two nodes (C, Nq) and (C", Nc) are connected via an edge if and only if 
C G C (or C C C) and there exists no other component C" with C C C" C C 
{C C C" G C). In the first case, we say that {C,Nc) is a child of {C'^Nc), 
whereas, in the second case, {C,Nc) is a parent of {C,Nc)- For brevity, we 

'^^More precisely, it is linear in \2.(2B)\ up to a factor that is polynomially bounded in logn, 
loglog(u)(B)), and loglog(crf’(2B)“^). If 2B contains no root, then there is only one iteration. 
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usually omit the integer Nc, and just refer to C as the nodes of T. Notice that, 
according to Theorem 4, the so obtained graph is indeed a tree rooted at B. A 
node C is called terminal if and only if it has no children. We further use the 
following definition to refer to some special nodes: 

Definition 3. A node (C, Nc) G T is called special, if one of the following 
conditions is fulfilled: 

• The node (C, Nq) is terminal. 

• The node {C,Nc) is the root ofT, that is, {C,Nc) = {B,4:). 

• The node {C,Nc) is the last node for which Bisection is called in the 
preprocessing phase of the algorithm. We call this node the base of T. 
Notice that the first part of the tree consists of a unique path connecting 
the root and the base of the tree. 

• For each child D of C, it holds that Z{D^) Z{C~^). 

Roughly speaking, except for the root and the base of T, special nodes either 
isolate a root of F or they are split into two or more disjoint clusters each 
containing roots of F. More precisely, from Theorem 4, we conclude that, for any 
two distinct nodes C,D G T, the enlarged regions Z{C~^) and Z{D'^) are either 
disjoint or one of the nodes is an ancestor of the other one. In the latter case, we 
have C'^ C D+ or Z3+ C C'*'. Since, for any two children Di and D 2 of a node C, 
the enlarged regions and are disjoint, we have X]?=i ^i^t) — 
where Di to are the children of C. Hence, since each Df contains at least 
one root, the fourth condition in Definition 3 is violated if and only if C has 
exactly one child D and Z{C'^) = Z(D+). The number of special nodes is at 
most 2 • (1 + \Z{2B)\) as there is one root and one base, at most \Z{2B)\ terminal 
nodes C with C B, and each occurrence of a special node, which fulfills the 
fourth condition, yields a reduction of the non-negative number '^c{\Z{C^)\ — 1) 
by at least one. The subdivision tree T now decomposes into special nodes 
and sequences of non-special nodes Ci,... ,Cs, with Ci D (72 D • • • D Cg, that 
connect two consecutive special nodes. The remainder of this section is dedicated 
to the proof that the length s of such a sequence is bounded by some value Smax 
of size 


Smax = (7 (logn-f loglog(w(H)-f loglog(CTF(2H) ^)) (12) 

= O (log (n ■ log{w{B)) ■ l6g(crF(2H)"^))) . 

For the proof, we need the following lemma, which provides sufficient condi¬ 
tions for the success of the NewtonTest. 

Lemma 6 (Success of NewtonTest). Let C = {Bi, ..., Bg^} be a non¬ 
terminal component with B \ C ^ let Be be the corresponding enclosing 
square of width w{C) and center m = me, and let A Ac = A(m,r), with 
r := jw{C), be the corresponding enclosing disk. Let zx,...,Zk be the roots 
contained in the enlarged component C'^, and suppose that all these roots are 
contained in a disk A" := A(m",r") of radius r” = In addition, 

assume that the disk A(m, 2^^°s""''^°Acr) contains none of the roots Zk+i,... ,Zn. 
Then, the algorithm CIsolate performs a quadratic step, that is, C is replaced 
by a single component C of width w{C') < . 
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Proof. We first argue by contradiction that 4A does not intersect any other 
component C , which implies that the first condition of the and in the if-clause 
in line 9 is fulfilled. If 4A intersects C, then the distance between C and C 
is at most 8r, and thus 2^^' < 8r as the distance between C and C" is at least 
max(2^'^, Hence, we conclude that the disk A(m,64r) completely contains 

2B' for some square B' of C". Since 2B' also contains at least one root and 
since each such root must be distinct from any of the roots zi,... ,Zk, we get a 
contradiction. 

According to our assumptions, each of the two disks A and 8A contains 
the roots zi,...,Zk but no other root of F. Hence, according to Lemma 4, 
T*(2A) = T*(4A) = k holds. Since we assumed C to be non-terminal, we must 
have k > 2, and thus the algorithm reaches line 11 and the NewtonTest is 
called. We assumed that C does not entirely cover the initial square B, hence, 
in a previous iteration, we must have discarded a square of width 2^^ or more 
whose boundary shares at least one point with the boundary of C. Hence, we 
can choose a point in such a square as the point xc G B\C in the NewtonTest 
such that the distance from xc to C is equal to 2^^~^ and such that the distance 
from Xc to the boundary of B is at least 2^'^“^. Notice that also the distance 
from Xc to any other component C’ is at least 2^'^“^, and thus the distance from 
Xc to any root of F is at least 2^'^“^, which is larger than or equal to as C 
consists of at most 9n squares. From our assumptions, we thus conclude that 

T 

\xc — m"| < 4r and \xc — Zi\> —— for i < k, 

27n 

and 

\xc — Zi\> 2^'^n^Nc • r — 4r > 2^®n^A'c ■ r for i > k. 

Using the fact that ^ ^i^h F{x) 0, we can bound 

the distance from the Newton iterate x'q as defined in (11) to the “center” m" of 
the cluster of roots: 


1 {xc - m")F'{xc) 
k F{xc) 



Xc — rn" 
xc - Zi 



- 1 


< 


1 

k 


E 


Zi — m" 
xc - Zi 
n — 


rl{27n) k 


+ E 

i^k 


Xc — rn" 

xc - Zi 


k 4r 
^ rF2^'^Ncr 


^ 1 -A \zi - m"\ ^ \xc - m"\ 
-k^^\xc-z,\ ^^\xc-z,\ 

27nr Far 1 

ra?2‘^°Nc ^ rn^2^^Nc ~ 2^‘^nNc' 


Hence, there is an e £ C, with |e| < 2 wbvF> k 'f{xc) = 1 + e. 

This implies that and thus the NewtonTest must 

reach line 2 as Algorithm 3 must return True or Undecided. With x'q = 
xc — k ■ it further follows that 


\m" — x'q\ = \m" — Xc\- 1 — 
e(m" - xc) 


I {xc-m.")F'{xc) 
k F(xc) 


1 -I- £ 


< 


4r 


< 


= \m - xc\ 
2^c 

< 


1 - 


1 + £ 


2i3nAc “ 2iinAc 128Ac ' 
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We can therefore bound 


oic 

\i'c - < \x'c - x'c\ + \x^G - m''\ < + Wc - m"\ 

2ic 2'c 2^^ 

- 64:Nc ^ 128Nc ^ 82Nc' 

Since the distance from m" to any of the roots zi,..., is also smaller than 
r" < ^2Nc ’ conclude that the disk A(a;p, contains all roots zi,... ,Zk- 

Hence, we conclude that A' := A(a;p, §]^) is ( 5 , |)-isolating for the roots 
zi,... ,Zk, and thus T*(A') = k must hold according to Lemma 4. This shows 
that we reach line 11 and that the NewtonTest returns SUCCESS. □ 


In essence, the above lemma states that, in case of a well separated cluster of 
roots contained in some component C, our algorithm performs a quadratic step. 
That is, it replaces the component C by a component C of width w{C') < < 

, which contains all roots that are contained in C. Now, suppose that there 
exists a sequence C'l,... ,Cs of non-special nodes, with Ci D • • • D Cg, such 
that Cs has much smaller width than Ci. Then, Ci contains a cluster of nearby 
roots but no other root of F. We will see that, from a considerably small (i.e., 
comparable to the bound in ( 12 )) index on, this cluster is also well separated from 
the remaining roots (with respect to the size of Ci) such that the requirements 
in the above lemma are fulfilled. As a consequence, only a small number of steps 
from Ci to Ci+i are linear, which in turn implies that the whole sequence has 
small length. For the proof, we need to consider a sequence {si)i = {xi,ni)i, 
which we define in a rather abstract way. The rationale behind our choice for 
Si is that, for all except a small number of indices and a suitable choice for 
Si, the sequence {si)i behaves similarly to the sequence (2^"^^, log log Aq-Ji. We 
remark that {si)i has already been introduced in [45], where it serves as a crucial 
ingredient for the analysis of the real root isolation method ANewDSC. 

Lemma 7 ([45], Lemma 25). Let w, w' € K“'" he two positive reals with w > w', 
and let m € N>i be a positive integer. We recursively define the sequence 
(si)iGN>i := ((a:i,?T-i))iGN>i as follows: Let si = {xi,ni) := {w,m), and 




{ci ■ Xi, + 1) with an a G [0, if jf: > w' 

{Si ■ Xi,max{l,ni — 1 )) with a Si G [ 0 , ^j, if ffi < w', 


where Ni := 2^ ’ and i >1. Then, the smallest index io with Xig < w' is bounded 
by 8{ni -I- loglogmax(4, ^)). 

We are now ready to prove the claimed bound on the maximal length of a 
sequence of non-special nodes: 

Lemma 8. Let V = {Ci,Ni ),..., (Cg, Ng), with Ci D • • • D Cg, be a sequence 
of consecutive non-special nodes. Then, we have s < Smax with an Smax of size 


Smax = O (logn -k loglog(w(H)) -f log log(crF(B+) ^)) 
= O (log (n • log{w{B)) ■ \og{aF{2B)~^))). 
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Proof. We distinguish two cases. We first consider the special case, where V is 
an arbitrary sub-sequence of the unique initial sequence from the child of the 
root of the tree to the parent of the base of the tree; if there exists no non-special 
root in between the root and the base of the tree, there is nothing to prove. Due 
to Theorem 4, part (e), Cs consists of at most 9 • \Z{2B)\ squares. It follows 
that 2^® < 9n as Ci consists of at least 2^* squares. This yields s = O(logn). 

We now come to the case, where we can assume that each Ci is a successor of 
the base of the tree. In particular, we have B\Ci ^ W.l.o.g., we may further 

assume that zi,..., are the roots contained in the enlarged component . 
Since all Ci are assumed to be non-special, each Cf contains zi to z^ but no 
other root of F. Let Wi := w{Ci) be the width of the component Ci, r* := | • Wi 
be the radius of the enclosing disk := Ac^, and 2^* := 2^‘^i be the width of 
each of the squares into which Ci decomposes. Notice that, for each index i, the 
enlarged component C^ is contained in the disk 2 Ai of radius | • Wi, and thus the 
disk 2As of radius | • Ws contains the roots zi to z^. We now split the sequence 
V into three (possibly empty) subsequences Vi = {Ci,Ni ),..., (Ci,, Wi), P 2 = 
(C,,+i,W,+i),...,(C, 2 ,WJ, and V 3 = (C^^+i, W.+i),..., (C^, fV^), where ii 
and 12 are defined as follows: 

• ii is the first index with 2^^ > • Ai, • 2^'i. If there exists no such 

index, we set ii := s. Further notice that, for any index i larger than ii, 
we also have 2^^ > 2^ • Ni ■ 2^*, which follows from induction and 

the fact that 2^* and Ni are replaced by and Nf in a quadratic step. 

• 12 is the first index larger than or equal to ii such that the step from 12 to 
12 + 1 is quadratic and 2^” • 2^• Ni^ > 2^*2. If there exists no such 
index 12 , we set 12 := s. 

From the definition of 12 , it is easy to see that Vs has length bounded by O(logn). 
Namely, if 12 = s, there is nothing to prove, hence we may assume that the step 
from ^2 to ia + l is quadratic and 2^= > A^^- 2^-2 = 2-3log"-3i . 2 ^. 2 +!. 

Hence, we conclude that s — (*2 + 1) A 31ogn -|- 31 as ii is reduced by at least 1 
in each step. 

Let us now consider an arbitrary index i from the sequence V 2 - The distance 
from an arbitrary point in C^ to the boundary of C^ is at least 2 ^^~^ > 
23i°g"+3i . A, • 2^- > 22 '°s"+20 ■ Ni-Ti, where the latter inequality follows 
from Ti = ^Wi < I • 9n • 2^^ Since C^ contains only the roots zi,..., z^, this 
implies that the distance from an arbitrary point in Cf to an arbitrary root 
Zfc+i, ■ ■ ■ ,Zn is larger than 2^*°s”+^° ■ Ni ■ Vi. Hence, the second requirement 
from Lemma 6 is fulfilled for each component Ci with i > ii- Now, suppose 
that 2 ^’ ■ 23i°g"+32 . < 2 ^% then the roots zi to Zk are contained in a disk 

of radius | • 9n • 2^“ < 2-23-*°g” • A”^ . n, and thus also the first requirement 
from Lemma 6 is fulfilled. Hence, from the definition of 12 , we conclude that 
the algorithm performs a quadratic step if and only if 2^® • 2^ "+32 . ^ 2 ^'. 

We now define the sequence Si := (2^% log log-^ 1)1 where i runs from H to the 
first index, denoted for which 2^’'i < 2^’ ■ 2-3iog™-32^ Then, according to 
Lemma 7, it holds that — ii < 8 {m + loglogmax(4, ^)), with w := 2^h , 
m := log log Ai,, and w' := 2^“ ■ 23iog"+32^ Theorem 4 (g) yields that m = 
0(logldg(w(H)) -I- logldg(CT_F(2H)“^). Hence, since 12 — i'l < 31ogn -I- 32, we 
conclude that h — ii < 0(\ogn + loglog(w(H)) -I- loglog(cr_F(2H)“^)). 
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It remains to show that the latter bound also applies to ii. From the 
upper bound on the numbers Ni, it follows the existence of an Wmax of size 
0(logn + logl6g(ri;(;B)) + loglog(CTF(2yB)“^)) such that each sequence of con¬ 
secutive quadratic steps has length less than Wmax) and such that after Wmax 
consecutive linear steps, the number Ni drops to 4. Since the number ii de¬ 
creases by at least 1 in each step, there exists an index i' of size O(logn) 
such that • 23'°s"+34 < 2 ^ 1 . Now, if the sequence C^, C^+i,... starts with 
Wmax or more consecutive linear steps, we must have = 4, and thus 

2^i'+™max . Hence, we conclude that ii < / -I- m-max in 

this case. Otherwise, there must exist an index i", with i' < i" < i' mmax) 
such that the step from i" to i” + 1 is quadratic, whereas the step from i" -f 2 is 
linear. Then, it holds that 


iVF '+2 = v'^i"+i = and 2^-"+2 < 2^-"+i 



< 2 


— 3 log n—32 


2^1 


which implies that ii < i" + 2 < i' + rumax + 1 = 0(logn -I- loglog(?c(;B)) -f 
loglog(crF(2S)“^)). Hence, the claimed bound on ii follows. □ 

We can now state the first main result of this section, which immediately 
follows from the above bound on Smax and the fact that there exists at most 
2 • {\Z{2B)\ + 1) special nodes: 

Theorem 5. The subdivision tree T induced by CIsolate has size 


|r|<2-(|Z(2,B)| + l)-s,„ax (13) 

= O {\Z{2B)\ ■ log (n • \^{w{B)) ■ \^{aF{2B)~^))) . 

If B contains all complex roots of F, and if \og{w{B)) = 0(rF -l-logn)/'^ then 
the above bound writes as 


O (n • log (n • Ff • log((T^^))) . (14) 

We can also give simpler bounds for the special case, where our input 
polynomial has integer coefficients. Suppose that f{x) G Z[x] has integer 
coefficients of bit size less than t. We first divide / by its leading coefficient 
lcf(/) to obtain the polynomial F := f / lcf(/), which meets our requirement from 
(1) on the leading coefficient. Then, we have Ff = 0(r) and ap = 
e.g. see [62] for a proof of the latter bound. Hence, we obtain the following 
result: 


Corollary 3. Let f be a polynomial of degree n with integer coefficients of bit 
size less than t, let F := //lcf(/), and let B be a square of width ^ 

Then, the algorithm CIsolate (with input F and B) uses 

0{\Z{2B)\ ■ log(nT)) = 0{n ■ log(nT)) 
iterations to isolate all roots of F that are contained in B. 

The above results show that our algorithm performs near-optimal with respect 
to the number of components that are produced by the algorithm. In addition, 

^^Notice that we can compute such a square B with 0{n^Tp) bit operations; see Section 2. 
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since each component consists of at most 9n squares, we immediately obtain an 
upper bound for the total number of squares produced by the algorithm that 
exceeds the bound from (14) by a factor of n. Indeed, we will see that the actual 
number of squares is considerably smaller, that is, of size 0(\Z(2B)\ ■ Smax • logn), 
which exceeds the bound in (14) only by a factor logn. For the proof, we consider 
two mappings cj) and if}, where (j) maps a component C = {Bi, ..., i?sc} fo ^ root 
Zi G and ip maps a square Bj to a root Zi G 4Bj n C’*'. The claimed bound 
for the total number of squares then follows from the fact that we can define 
(p and Ip in a way such that the pre-image of an arbitrary root Zi G 2B (under 
each of the two mappings) has size 0(smax • logn). The rest of this section is 
dedicated to the definitions of p and ip and the proof of the latter claim. In what 
follows, we may assume that 2B contains at least one root as, otherwise, all four 
sub-squares of B are already discarded in the first iteration of the preprocessing 
phase. 

Definition 4. For a root ^ G 2B, we define the canonical path of ^ as the 
unique path in the subdivision tree 7b that consists of all nodes C with ^ G C'^ . 

Notice that the canonical path is well-defined as, for any two nodes Ci and 
C 2 , either and are disjoint or one of the two components contains the 
other one. We can now define the maps p and ip: 

Definition 5 (Maps p, ip). Let C = {Bi ,..., Bsc} ® node in the subdivision 
tree Tb, and let B := Bj be an arbitrary square in C. Then, we define maps p 
and Ip as follows: 

(p) Starting at C, we descend in the subdivision tree as follows: If the current 
node D is a non-terminal special node, we go to the child E that minimizes 
\Z{E'^)\. If D is terminal, we stop. If D is non-special, then there is a 
unique child of D to proceed with. Proceeding this way, the number \Z{D~^)\ 
is at least halved in each non-terminal special node D, except for the base 
node. Hence, since any sequence of consecutive non-special nodes has length 
at most Smax; it follows that after at most Smax • (logr(l'2'(C'“'')|)] -b 1) < 
Smax ■ (logn -b 2) many steps we reach a terminal node F. We define p{C) 
to be an arbitrary root contained in Z{F'^). 

(ip) According to part (d) of Theorem f, the enlarged square 2B contains at least 
one root Now, consider the unique maximal subpath P( = Ci, C 2 ,..., Cs 
of the canonical path P^ that starts at Ci := C. If s < |"log(18n)], we 
define ip{B) := Otherwise, consider the component C := C'pog(i8n)] 
and define ip{B) := p{C'). 

It is clear from the above definition that p{C) is contained in C'"*' as each 
root contained in the enlarged component F'^ corresponding to the terminal 
node F is also contained in . It remains to show that ip{B) G 4:B n C+. If 
the length of the sub-path P( is |"log(18n)] or less, then ip{B) = f, G 2B, hence, 
there is nothing to prove. Otherwise, the squares in C have width less than 
™ . Since C can contain at most 9n squares, we conclude that w{C') < , 

and since ^ is contained in B+ as well as in [C')^ , we conclude that (C")+ C 4B, 
and thus p{B) = P{C') G 4B O {C')+ C 4B n C+. 

Now, consider the canonical path P^ = C'i,...,Cs, with Ci := B, of an 
arbitrary root ^ G 2B. Then, a component C can only map to ^ via p if 
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C = Ci for some i with s — i < Smax • (log \Z{2B)\ + 1). Hence, the pre- image 
of ^ has size 0(smax • log |-Z(2H)|). For the map ■(/>, notice that a square B 
can only map to ^ if H is contained in a component C = Ci for some i with 
s — i = Smax • (log|-Z(2H)| + 1) + |’log(18n)]. Since, for each component Ci, 
there exist at most a constant number of squares B' G Ci with ^ G 4B', we 
conclude that the pre-image of ^ under ijj is also of size 0(smax • ldg2(2H)). 
Hence, the total number of squares produced by our algorithm is bounded by 
0{\Z{2B)\ ■ Smax • log \Z{2B)\). We summarize: 

Theorem 6. Let ^ G 2B be a root of F eontained in the enlarged square 2B. 
Then, with mappings cj) and ip as defined in Definition 5, the pre-image of f, 
under each of the two mappings has size 0(smax ’ log \Z{2B)\). The total number 
of squares produced by the algorithm CIsolate is bounded by 

0(smax • \Z{2B)\■ log \Z{2B)\) = d{n- log {log{w{B)) ■ \og{aF{2B)-^))) 

We can also state a corresponding result for polynomials with integer coeffi¬ 
cients: 

Corollary 4. Let f G Z[x] and F := //lcf(/) he polynomials and B be a square 
as in Corollary 3. Then, for isolating all roots of f contained in B, the algorithm 
CIsolate (with input F and B) produces a number of squares bounded by 

0{\Z{2B)\ ■ log(nr) • logn) = 0(nlog^(nr)). 

5.2 Bit Complexity 

For our analysis, we need to introduce the notion of a square or component being 
weakly centered and centered: 

Definition 6. We say that a square B of width w := w{B) is weakly centered, 
if min{| 2 ;| : z G B} < A ■ w. Similarly, we say that a square is centered if 
min{|z| : z G B} < w/A. In addition, we define a (weakly) centered component 
to be a component that contains a (weakly) centered square. 

Notice that the child of a component that is not weakly centered can never 
become weakly centered. Hence, it follows that the set of weakly centered 
components forms a subtree Twcent of the subdivision tree T that is either empty 
or contains the root component B; see Figure 4 for an illustration. 

Moreover, let C and C be siblings in Twcent and let w and w' be the 
sizes of the boxes in the components C and C, respectively. We already 
argued that the distance between C and C is at least max{r(;, w'}. W.l.o.g. let 
min{|z| : z G C} < min{|z| : z G C'}, then the distance of C' to the origin is at 
least w'/2, and thus C' is not centered. We further conclude that a descendant 
of depth 3 of C is not weakly centered because the width of the boxes in the 
component is at least halved in each step. It follows that each path in Twcent 
consisting of only weakly centered components has length at most 3. From this 
observation, we conclude that the subtree Twcent has a very special structure. 
Namely, it consists of only one (possibly empty) central path P = Ci,... ,Ci 
of all centered components, to which some trees of depth at most 4 of weakly 
centered components are attached, see Figure 4. Since there can only be a 
constant number of disjoint not weakly centered squares of the same size, it 
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further holds that the degree of each node is bounded by a constant. Hence, 
each of the attached trees has constant size, and each node in the tree contains 
at most a constant number of weakly centered squares. 

Notice that not weakly centered components C € 'T\ T^cent have the crucial 
property that any two points in C'^ have absolute values of comparable size; 
see Lemma 9. This is not true in general for (weakly) centered components as 
the size of C might be very large whereas the distance from to the origin is 
small. 

Lemma 9. If B is a not weakly centered square, it holds that 

maxldg( 2 ) < 5 + min Idgz. 

zG4B 2G4B 

Moreover, if C is a not weakly centered component, then it holds that 
max Idg(z) < log(64n) + min log 2 . 

zGC+ zGC+ 

Proof. We first prove the claim for a not weakly centered square B. Let z := 
argmin{|z| : z G 45} and z := argmax{|z| : z G 45}. By definition the distance 
of 5 to the origin is at least 4w, where w denotes the size of 5. Moreover, with 
X := argmin{|z| ; z G 5}, we get |z| > |a:| — |a; — z| > iw — Vl2.5ic > w/A. Thus 

1^1 ~ \z\ < 1^ — z| < Ay/2w < 16\/2 ■ |z|, 
and the statement follows. 

It remains to prove the claim for a not weakly centered component. Let 
z := argmin{|z| : z G C+}, z := argmax{|z| : z G C+}, and let 5 C C be a 
square in C such that z G 45. Then, as above, it follows that z > w/4 and since 
C contains at most 9n squares it holds that 

1^1 ~ \z\ < |z — z| < ■ {9nw + w) < \f2 ■ lOnw < 57n|z|, 

and thus max 2 gc+ Idg(z) < log(64n) + minj,gG?+ Ibgz. □ 

In the previous section, we introduced mappings (j) and that map com¬ 
ponents C and squares 5 to roots contained in C'*' and 45 n C'*', respectively, 
such that the preimage (under each of the two mappings) of each root has size 
at most 0(si„ax ■ log |Z(25)|) = 0(smax • logn), with 

Smax = 0(log(nldg(w(5))log((Ti2(25)“^))) 

as defined in Lemma 8. The crucial idea in our analysis is to bound the cost for 
processing a certain component C (square 5) in terms of values that depend only 
on the root (j){Ct) {ijj{B)), such as its absolute value, its separation, or the absolute 
value of the derivative F' at the root; see Lemma 13. Following this approach, 
each root in 25 is “charged” only a small (i.e. logarithmic in the “common” 
parameters) number of times, and thus we can profit from amortization when 
summing the cost over all components (squares). For each not weakly centered 
component C, the width of C and the absolute value of any point in C is upper 
bounded by 64n • \4>{C)\, which allows us to bound each occurring term ldg4(;(C') 
by 0{\ogn -|- log |(()(C)|). However, for weakly centered components (squares), 
this does not hold in general, and thus some extra treatment is required. For 
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Figure 4: A possible subdivision tree of CIsolate on an arbitrary input box. 
Red nodes (squares and diamonds) correspond to centered components, blue 
nodes (large circles) to weakly centered components and black nodes (small 
circles) to not weakly centered components. The subtree Tiwcent consists of all red 
and blue nodes. Note that 7(vcent consists of one path of centered components 
(marked by the red edges) and trees of depth at most 4 attached to it (consisting 
of blue nodes). The rectangular red nodes correspond to those components 
Cij,..., Ci^ on the central path that are split nodes according to Definition 3, 
that is, components for which Z(C(*i) D The difference of these two 

sets contains the roots that we map the (weakly) centered components to. For 
instance, the centered box C' is mapped to a root (j){C) that is contained in 
Z(C')^^j^) \ Z(C)^). All weakly centered boxes that have C as their first centered 
predecessor, as for example C, are mapped to the same root 


this, we will introduce slightly modified mappings </> and ip that coincide with 
(p and Ip on all not weakly centered components and squares, respectively, but 
map a weakly centered component (square) to a root of absolute value that is 
comparable to the size of the component (square). In the next step, we will show 
that this can be done in a way such that the pre-image of each root is still of 
logarithmic size. We give details: 

Let ii, ..., is, with 1 < ii < ^2 < .. . < is < i—1, be the indices of components 
in the central path P = Ci, ... ,Ci such that Z{C^.) A Z(C^.^i) for j = 1,..., s. 
That is, Ci^ are the weakly centered components that are also special according to 
Definition 3. In addition, we say that Zq := max{|a;| : x G 2B} is the pseudo-root 
of2B, and define Z~^{2B) := Z{2B) U {zq} as the set consisting of all (pseudo-) 
roots in 2B. 

Definition 7 (Maps (p and ip). Let C be a component and B C C be a square 
contained in C. 

1 . If the component C is not weakly centered, we define <p{C) := <p{C). 

2. If the component C is weakly centered, let Ci G P be the first centered 
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predecessor of C in Twcent- If there exists no such Ci or if i G we 

define 4>{C) = zq. For i G (* 2 )^]; j be maximal with ij < i. Then, there 
exists a root ^ G Z(C^) \ We define 4>{C) := f. 

3. If B is weakly centered, we define ip{B) := 4>{C). Otherwise, we define 
^(B) := if{B). 

We derive the first crucial property of the mappings (j) and if: 

Lemma 10. It holds that 

max log z < log(64n) + log <?!>(C') and max log z < log(64n) + log if{B) (15) 

zGC+ zG4B 

for all components C (squares B). 

Proof. For a not weakly centered component C (square B C C'), this follows 
directly from Lemma 9 and the fact that (j){C) = 4>{C) {if{B) = if{B)) and that 
(fiC) G C+ {if{B) G 4BnC+). 

For a centered component C, we either have <f>{C) = zq or (j){C) ^ C'"*' due 
to the definition of (f. The first case is trivial, hence, we may assume that 
(t>{C) ^ C'^. Since C is centered, it contains a centered square B, and thus the 
distance of B to the origin is at most w/A. It follows that contains the 
disk of radius w/A around the origin, hence | 9 !>(C')| > w/4. Since the distance 
between any two points in C~^ is upper bounded by IQ'/^nw, we conclude that 
|z| < w/A + 10^/2nw < (1 + AOy/2) ■ \${C)\ < 64n • \^{C)\ for all z S (7+. The 
same argument further shows that |z| < 64n • |'0(^)1 for centered squares and 
all z e 4i?. 

It remains to show the claim for a weakly centered component C (square B) 
that is not centered. In this case, we either have '0(C') = zq or if{C) = 0(C"), 
where C is a centered component on the central path that contains C; see 
Definition 7. In the first case, there is nothing to prove. In the second case, 
we have have already shown that Inequality (15) holds for C', hence, it must 
hold for C as well. The same argument also applies to squares B in C that are 
weakly centered but not centered. □ 

Notice that, for a centered component C (square B), the image of the 
corresponding mapping 0 (if) may no longer be contained in the enlarged 
component C~^ (enlarged square AB), as it is the case for the mappings (f and if. 
However, it still holds that the preimage of its (pseudo-) root under each of the 
two mappings cf and if is of small size: 

Lemma 11. Let ^ G .Z~^(2B). Then, the preimage of under (f and if has size 
at most 0(sinax • lbg|Z(2;B)|) = 0(si„ax • logn). 

Proof. Since (f coincides with (f on all components C ^ T^cent, it suffices to show 
the claim for the restriction (f\T„cBat of 4’ to the components C G Pwcent that are 
weakly centered. Let ^ = if{C), with C G T^cent, be an arbitrary root contained 
in the image of let Ci be the first predecessor of C that is central 

and located on the central path. Then, there exists a j G {1, ... — 1} with 

ij < i < ij+i, and we have ^ G Z{C('.) \ Z(C)*l^j^). In addition, C is connected 
with Ci via a path of constant length as the distance from C to the central path 
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on Twcent is bounded by a constant and there cannot be more than 3 consecutive 
components that are weakly centered but not centered. Since there exists at 
most Smax Components on the central path between Ci- and Cij+i, it follows 
that the number of components C S Twcent that are mapped to ^ is bounded by 
O(Sinax)- The same argument applies to the special case, where C is mapped to 
the pseudo-root zq. Also, from the same argument and the definition of "0, it 
further follows that there can be at most O(sniax) many weakly centered squares 
that are mapped to the same (pseudo-) root, since each component contains at 
most constantly many weakly centered squares. □ 

We can now start with the bit complexity analysis of CIsolate. Let C = 
{Bi,..., Ssj,} be any component produced by the algorithm. When processing C, 
our algorithm calls the T»-test in up to three steps. More specifically, in line 9 of 
CIsolate the T*(2Ac) and the T*(4Ac)-test are called. In the NewtonTest, 
the T*(A')-test is called, with A' as defined in line 4 in NewtonTest. Finally, 
in Bisection, the T*(AB')“test is called for each of the 4 sub-squares B' into 
which each square Bi of C is decomposed. Our goal is to provide bounds for the 
cost of each of these calls. For this, we mainly use Lemma 5, which provides a 
bound on the cost for calling T*(A) that depends on the degree of F, the value 
Tp, the size of the radius and the center of A, and the maximal absolute value 
that F{x) takes on the disk A. Under the assumption that A has non-empty 
intersection with C, we may reformulate the latter value in terms of parameters 
(such as the absolute value, the separation, etc.) that depend on an arbitrary 
root contained in . 


Lemma 12. Let C be a component, and let A := A(to, r) be a disk that has non¬ 
empty intersection with C. If ZiC'^) > 1, then it holds that aF{zi) < n ■ 2^"^+® 
and 

max |F(z)| > 2“^®" • apizt) ■ \F'(zi)\ ■ {\Z{C'^)\ ■ maxi A)“^", 

A 

where Zi is an arbitrary root of F contained in C'^ and A := the ratio 
of the size of a square in C and the radius of A. If Z(C'^) = 0, then C 
is the component consisting of the single input square B, and it holds that 
max^gA |T'(-z)| > (2ui(S))“”“^. 


Proof. Notice that \Z{C^) \ = 0 is only possible if C = S as contains at least 
one root for each component C that is not equal to the single input square B. 
Hence, each point z G A n C has distance at least w{B)/2 to each of the roots 
of F, and thus |F(z)| > (2w(B}}-^-^. 

In what follows, we now assume that C~^ contains at least one root f. Let 
us first bound the separation of If there exists another root G C~^, then 
we must have api^ < |C - ?'l < (9 • \Z{C+)\ -h 1) • f • < „ . 2 ^c+^, where 

we used that each component C consists of at most 9 • \Z{C'^)\ squares, each 
of size 2^^. Now, let \Z{C^)\ = 1, and let C be the direct ancestor of C. 
When processing C, the NewtonTest failed as its success would imply that 
k := \Z{{C')'^)\ = \Z{C'^)\ > 2. In addition, since C is non-terminal, the disk 
8Ac/ must contain at least two roots as otherwise T*(2Ac/) as well as T*(4Ac/) 
would return 1. Hence, we have (Tp{f) < n • 2^c'+5 = ri ■ 2^'^+®. We conclude 
that, in any case, any root ^ G C~^ has separation ap{f,) < n ■ 2^'^+®. 

In the next step, we show that there exists an m! G AnC’*' whose distance to 
C is at most 2^'^“^ and whose distance to any root of F is at least • 
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Namely, due to our assumption, there exists a point p G AO C. Then, the two 
disks A(p, and A share an area of size larger than min(2^^“^, r)^, and 

thus there must exist an m' G A{p, 2^'^“^) n A C C~^ whose distance to any root 

of F is lower bounded by ^ 

Now, let Zi G C~^ be an arbitrary but fixed root of F. If Zj is a root not 
contained in (7+, then \m' — Zj\ > 2^"^“^, and 


\Zi - Zj\ 

\m' - Zj\ 


< 


\zi — m'\ + \m' — Zj 


< 1 + 


< 1 


m' — Zi 


TO' — Zj 


9 • n • 2^c+i 

2<^c-2 


= 1 + 72 • \Z{C+)\ < 2^ 


' n. 


If Zj is a root in C"*", then 

\zi — Zj\ 9712 ^*^+^ 
\m' — Zj\ ~ S 

Hence, we get 


18 • 71^/2.2^c+i 

min(2^c-2^ f'j 


< 2® • 77®/^ • maxi A. 




|+(m')| = \Fn\ • |to' - Zj\ = \F'{zi)\■ \m' - +| • H 1- rV 

> \F'izi)\ ■ \m' - z,| • (2®77®/2 maxi 

> \F {Zi)\ ■ --(2®77®/^ maxi A) 

2y/n 

'' A ^ P'’”’'""''"' 

> 2“^®” • apizi) • \F'{zi)\ • (77 • maxi A)”^"^ 


where, in the second to last inequality, we used that apizi) < 2^^+® • n. □ 

We can now bound the cost for processing a component C: 

Lemma 13. When processing a component C = {Hi,..., Bs^} with |Z(C'+)| > 
1, the cost for all steps outside the NewtonTest are bounded by 

0(77 • {n\og(j){C) + \og{ap{(j){C))~'^) + \og{F'{(j){C))~'^))) + 

SC 

d{n- (Y2 n log if {Bi) + Tp + \ogap{il;{Bi))~'^ + lbgF'(V'(H,))"^)) (16) 

bit operations. The cost for the NewtonTest is bounded by 

0(77(77log UC) + log ap{f{C))-^ + logH'(0(O))-i + |Z(0+)| • log Nc)) (17) 

bit operations. If O"*" contains no root, then C is the component consisting of 
the single input square B, and the total cost for processing C is bounded by 
0{n{Tp + n\og(w{B),w{B)~^))) bit operations. 

Proof. We start with the special case, where O^ contains no root. Notice that 
this is only possible if O = H and 2B contains no root. Hence, in this case, 
the algorithm performs four T*(A)-tests in the preprocessing phase and then 
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discards B. Due to Lemma 5 and Lemma 12, the cost for each of these tests 
is bounded by 0{n{Tp + n\og{w{B),w{B)~^)) bit operations. Hence, in what 
follows, we may assume that C'"*' contains at least one root. We first estimate 
the cost for calling T* on a disk A = A{m,r), where A = 2Ac, A = 4Ac, or 
A = Api, where B' is one of the four sub-squares into which a square Bi is 
decomposed. If A = 2Ac or A = 4Ac, then lc)g(TO,r) = Oilogn + ldg(^(C')), 
and thus the cost for the corresponding tests is bounded by 

0(n • (nldg(^(C') + Tp+ F'{(j){C))~^)) 

bit operations, where we again use Lemma 5 and Lemma 12. For A = A^/, we 
have log(TO,r) = 0(\ogn + \og^p{C)), and thus Lemma 5 and Lemma 12 yields 
the bound 

0(n • (nl^ipiB,) + Tp + log api-ipiBi))-^ +\^F'{'tp{Bi))~^)) 

for processing each of the four sub-squares B' C Bi into which Bi is decomposed. 
Hence, the bound in (16) follows. We now consider the NewtonTest: In line 2 
of the NewtonTest, we have to compute an approximation x'q of the Newton 
iterate x'p, such that \x'q — x'(~.\ < For this, we choose a point xc G B\C 

in line 12 of CIsolate whose distance to C is 2^'^“^ and whose distance to 
the boundary of B is at least 2^'^“^. Since the union of all components covers 
all roots of F that are contained in B, and since the distance from C to any 
other component is at least 2^*^, it follows that the distance from xc to any 
root of F is larger than 2^'^“^. With A := A(a;c, 2^'^“^), it thus follows that 
\F{xc)\ > 2“” • max^gA |.P(z)|j and using Lemma 12, we conclude that 

IdgF(xc)"^ = 0{n\ogn F\^ap{(j){C))~^ +logF'{(j){C))~'^). (18) 

It follows that the cost for calling Algorithm 3 in Line 1 of the NewtonTest is 
bounded by 

0(n • (ldgF(xc)“^ -I- Tp -I- nldgxc)) 

= d{n ■ {tp + log a p{(j){C))~^ + l^ F'{(j){C))~^ +nlog^{C))) 

bit operations, Namely, Algorithm 3 succeeds with an absolute precision bounded 
by 0(ldg |T(xc)|~^) and, within Algorithm 3, we need to approximately evaluate 
F and F' at the point xc to such a precision; see also [20, Lemma 3] for the cost of 
evaluating a polynomial of degree n to a certain precision. If we pass line 1, then 
we must have |F'(xc)| > ^ gj:(c) ^ ■ Hence, in the for-loop of the NewtonTest, 
we succeed for an L of size 0(ldg F{xc)~^ + n log w)!!!) — t'c + log Nc), and thus 
the cost for computing x'p. is bounded by 

0{n- (ldgCTF(<?!>(C'))”^ + logF\(j){C))~^ + log^{C) + logNc)), 

where we again use (18) and the fact that ap{(j){C)) < n ■ 2^"^+® and ldgr(;(C') = 
0{logn -|- log(^(C')). It remains to bound the cost for calling T* on A' := 
A{x'c, I • in the NewtonTest. Again, we use Lemma 5 and Lemma 12 to 
derive an upper bound of size 

0{nlogn + l^ap{(j){C))~^ + logF'((()(C'))“^ -h |-Z(C'+)| • log Aq) 
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for log(max 2 gA' |^(-2)|) and thus a bit complexity bound of size 

d{n ■ (log^(C) + log+ log F\^{C))-^ + \Z{C+)\ ■ log Nc)) 

for T*(A'). This proves correctness of the bound in (17). We finally remark that 
the cost for all other (mainly combinatorial) steps are negligible. Namely, the 
bit size be of a square in a component C is bounded by 0(lbg(u>(C'), ^((7)“^) + 
logn) = 0{\ogn + log a F{(t>{C))~^ + lbg(/)(C')). Hence, each combinatorial step 
outside the NewtonTest, such as grouping together squares into maximal 
connected components in Line 8 of Bisection, needs 0{n) arithmetic operations 
with a precision 0{bc), and thus a number of bit operations bounded by (16). 

In the NewtonTest, we need to determine the squares Bi j of size /Nc 
that intersect the disk A'. This step requires only a constant number of 
additions and multiplication, each carried out with a precision bounded by 
0(log(w(C), + logn + log Ac). Hence, the cost for these steps is 

bounded by (17). □ 

When summing up the bound in (16) over all components C produced by 
the algorithm, we obtain the bit complexity bound 

d(^n- (n- (lZ(2B)l+logMeaF(2B)+log(w(B),w(B)-^))+TF ■ |Z(2S)| 

+ (l6gaF(zi)~^ 

z,eZ(2B) 

for all steps outside the NewtonTest. Here, we exploit the fact that the 
preimage of each (pseudo-) root in Z'^{2B) under each of the mappings <f>, </), tp, ip 
has size 0(sniax • logn), and that |zo| > w{B)/2 for the pseudo-root zq € Z^{2B). 
If we now sum up the bound (17) for the cost of the NewtonTest over all 
components, we obtain a comparable complexity bound; however, with an 
additional term n • "Y/c \^{^'^) \ ' log Ac, where the sum is only taken over the 
components for which the NewtonTest is called. The following considerations 
show that the latter sum is also dominated by the bound in (19). 

Lemma 14. Let 7new gT he the set of all components in the subdivision tree 
T for which the NewtonTest is called. Then, X^cgTnew ' log-l^C 

bounded by 

d{n-{\ogw{B) + \ogMeaF{2B) + \Z{2B)\)+TF-\Z{2B)\+ ^ lbgF'(z,)”^)- 

Zi&Z(2B) 

Proof We define := {C G Tnew : Ac = 4} and := {C G Tnew : 
Ac > 4}. Then, we have 


E \2{C^)\-^ogNc< E 2n< ^ 2n < 2n • \Z{2B)\ 


ceTz 


cer 


hence it remains to consider only the components C G 7^^^. For such a 
component, let anc*(C') gT he the last ancestor for which the NewtonTest 
succeeded. According to Theorem 4 (f), we have Ac < 4 • {w{anc*{C))/w{C)Y, 
and thus \Z{C^)\ ■ log Ac < 2 • |Z(C'+)| • (1 + log w(anc*(C')) — logw(C')). In 
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order to bound the later expression in terms of values that depend on the roots 
of F, let Zi G (7+ be an arbitrary root in . Then, assuming w{C) < 1, we 
obtain 


\F' {Zi)\ = \an\ ■ \Z^-ZJ\ \Zi-Zj\ 

j^i-.Zj&Z{2Ac) j:Zj^Z{2Ac) 

l^n I 

< 22 "+^^u>(C)I^( 2 Ac)|-i . inaxi(z,)”, 

where the last inequality follows from Landau’s inequality [56, Theorem 6.31], 
that is, Meap < ||p ||2 for any complex polynomial p G C[x], and \\F{zi — a ;)||2 < 
\\Fizi - x)||i < 2^^ maxi(zi)”2"+7 For a more detailed derivation of the latter 
see [45, Lemma 22]. Furthermore, using that Z{2Ac) contains at least two roots 
(as the NewtonTest is called) and that C'"*' C 2Ac, yields 

|F'(z,)l < 22’^+"^maxi(z,)"w(C)^^^ < 22"+"^ maxi(z,)”M^(C)^^. 

With Zi := (j){C) and lbg(/)(C') < log(64n) + Idg^(C'), we conclude that 

— |Z(C'''")| • logw(C') < 6nlog(64n) + 2tf + 2nlbg(^(C') + logF'{(j>{C))~^. 

(20) 

Notice that the latter inequality is trivially also fulfilled for a component C with 
w{C) < 1. In addition, it holds that 

|Z(C'''')| • logry(anc*(C')) < |Z(C^)| • (ldg(^(anc*(C)) + log(128n)). (21) 

The sum of the term at the right side of (20) over all C can be bounded by 

0(smax • (|-Z(2F)| • tf + ^ nlog{zi)+ ^ logF'(2;i)“^)) 

Zi^Z*{2B) ZieZ(2B) 

= 0(smax • (|-Z(2S) I • tf+ n log w(S) + n log MeaF(2S) + ^ log F'(z^)~^)), 

ZieZ(2B) 


as the preimage of each (pseudo-) root Zi G 2B (under (p and (j>) has size at 
most Smaxlogn. We may further omit the factor Smax in the above bound as 
logl6gtTF(2:i)“^ = 0(log(TF + nlog{zi) + logF'{zi)~^)) for an arbitrary root Zi 
of F, and thus 

Smax = 0(logn -b logldg'u;(F) -b loglbgCTF(2F)“^) 

= 0(logn-b logldg w(F)-b logTF-b logl6gMeaF(2F)-b log ^ lbgF'(zi)“^). 

ZieZ(2B) 

It remains to bound the sum of the term |Z(C'“'")| Tog (/)(anc*(C')) on the right 
side of (21) over all C G Let Tanc- := {C* G T : 3(7 G T with anc*((7) = 

C*}. Then, for a fixed C* G Tunc*, it holds that any (7 G with anc*((7) = 
C* is connected with C* in T via a path of length at most s := loglbgr(;(F) -b 
loglogtTF(2S)“^. Namely, we have already shown that loglogiVc < s for all 
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C G T, and thus Nq > 4 implies that the path connecting C and C* must have 
length at most s. Hence, it follows that 


|Z((C)+)| •logw(anc*(C)) 

C6r>,t:anc*(C)=C* 

^ |Z((C)+)|-log«;(C'*) 

Cer>,t:anc*(C)=C* 

< (ldg<^(C*) + log(128n)) ■ ^ \2{C+)\ 

CGr>/„:anc*(C)=C- 

< s • |Z((C'*)+)| • (l6g.^(C*) + log(128n)) 

< s • |Z(2H)| • (log<^{C*) + log(128n)), 

where we use the fact that, for any two components Ci,C 2 in the above sum, 
either n C 2 = 0, C or C C^. We conclude that 

l^(C'’^)l ■ logw(anc*(C')) 

<s-|Z(2H)|- ^ (log(^(C'*)+log(128n)) 

c*Gr,„c* 

= 0(smax • s • n log n • ^ (log(2j) + log(128n)) 

ZiGZ+(2B) 

= 0(smax ■ s-n - (\ogw{B) + logMea^(2H) + |Z(2H)|)) 

= 0(n • (l6gw;(H) + logMeaF(2H) + |Z(2H)| + log ^ logF'(zi)“^)). □ 

ZieZ{2B) 

Let us summarize our results: 

Theorem 7. Let F be a polynomial as defined in (1) and let B G C be an 
arbitrary axis-aligned square. Then, the algorithm CISOLATE with input B uses 

d{n ■ {n ■ {Z{2B) +\ogMea.F{2B) +\og{w{B),w{B)-^)) + tf ■ \Z{2B)\ (22) 

+ (McTFizi)-^+logF'{zi)-^))), 

Zi&Z(2B) 

bit operations. As input, the algorithm requires an L-bit approximation of F 
with 

6 {n ■ {Z{2B) + \^Mea.F{2B) + \og{w{B),w{B)-'^)) + tf ■ \Z{2B)\ (23) 

+ X! (log 0 -^( 2 :*)"^+lbgF'(zi)"^)). 

Zi^Z{2B) 

Proof. The bound (22) on the bit complexity follows immediately from Lemma 13, 
Lemma 14, and the remark following Lemma 13. The bound (23) on the precision 
demand follows directly from our considerations in the proof of Lemma 13 and 
Lemma 5. □ 

Notice that the above complexity bounds are directly related to the size of 
the input box B as well as to parameters that only depend on the roots located 
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in 2B. This makes our complexity bound adaptive in a very strong sense. In 
contrast, one might also be interested in (probably simpler) bounds when using 
our algorithm to isolate all complex roots of FI In this case, we may first 
compute an input box B of width w{B) = 2'"+^ that is centered at the origin. 
Here, T G N>i is an integer bound for Tp- with T = Tp- + O(logn). We have 
already argued in Section 2 that such a bound F can be computed using 0(n^r p) 
bit operations. Such a box B contains all complex roots of and thus running 
CISOLATE with input B yields corresponding isolating disks. Hence, we obtain 
the following result: 

Corollary 5. Let F be a square-free polynomial as in (1). Then, for isolating 
all complex roots of F, CISOLATE needs 


= d{n ■ {n^ + nlogMeaF + logT'(zi) ^)) (24) 

i=l 

= 0{n ■ {n^ + nlogMeaF + IbgDisc)^^)). (25) 

bit operations, where DIscf := Yli<i<j<n(^j ~ discriminant 

of F. As input, the algorithm requires an L-bit approximation of F with 

0(n^ + nlogMeaF + log Disc)^^) (26) 

Proof. The above bounds follow directly from the bounds in (22) and (23), and 
the fact that nlog{w{B),w{B)~^) nrp = 0{n'^ TnlogMeaF), 

n n 

^ldgcrF(zi)“^ = 0(n^ +nl6g(MeaF) +y^ldg F'{zi)~^), 

2 = 1 2=1 


and 

n 

\og F' = 0{r? + nlogMeaF + IdgDisc)^^); 

1=1 

e.g., see [45, Section 2.5] and the proof of [45, Theorem 31] for proofs of the 
latter bounds. □ 

Again, we provide simpler bounds for the special case, where the input 
polynomial / has integer coefficients: 

Corollary 6. Let / G Z[a:] be a square-free integer polynomial of degree n with 
integer coefficients of bit size less than t, let F := //lcf(/), and let B C C be 
an axis-aligned square with 2“'^*^'^^ < w{B) < 2'^*^'^^. Then, CIsolate with input 
B needs 0(nA + n^r) bit operations. The same bound also applies when using 
CIsolate to compute isolating disks for all roots of f. 

Proof. The claimed bound follows immediately from (25) and the fact that 
DiscF = lcf(/)^”“^ • Disc/ > For the second claim, we may simply 

run CIsolate on a box B of width 2"^+^ centered at the origin. According 
to Cauchy’s root bound, B contains all roots of /, and thus CIsolate yields 
corresponding isolating disks. □ 
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6 Conclusion 


We proposed a simple and efficient subdivision algorithm to isolate the com¬ 
plex roots of a polynomial with arbitrary complex coefficients. Our algorithm 
achieves complexity bounds that are comparable to the best known bounds 
for this problem, which are achieved by methods based on fast polynomial 
factorization [15, 28, 32]. Compared to these methods, our algorithm is quite 
simple and uses only fast algorithms for polynomial multiplication and Taylor 
shift computation but no other, more involved, asymptotically fast subroutines. 
Hence, also by providing a self-contained presentation and pseudo-code for all 
subroutines of the algorithm, we hope that there will soon be implementations 
of our method. So far, we have not discussed a series of questions concerning an 
efficient implementation, including heuristics and filtering techniques to speed 
up the computations in practice. In particular, we remark that Graeffe iteration 
is not well suited for implementations that are restricted to single- or double 
precision arithmetic. This is due to the fact that, after only a few iterations, 
it produces very large intermediate values with exponents in the correspond¬ 
ing floating point representations that are outside of the allowed range of the 
IEEE 754 specifications. However, when using multi-precision arithmetic, which 
would be a natural choice when implementing our solver, this issue is no longer 
relevant. In fact, preliminary tests on a few examples show that, when using 
multi-precision arithmetic, the relative error in the T^-test stays quite stable 
even after a number of Graeffe iterations.^® Hence, we are confident that a 
careful implementation of the T* will turn out to be efficient. 

Another possible direction of future research is to extend our current Newton- 
bisection technique and complexity analysis to the analytic roots algorithm 
in [61]. See ]52] for an alternative approach for the computation of the real roots 
of analytic functions obtained by composing polynomials and the functions log, 
exp, and arctan. 

At the end of Section 4.2, we sketched how to use our algorithm to isolate 
the roots of a not necessarily square-free polynomial for which the number 
of distinct complex roots is given as additional input. Also, we may use our 
algorithm to further refine the isolating disks for the roots of a polynomial 
in order to compute L-bit approximations of all roots. There exist dedicated 
methods ]20, 28, 32, 35, 36] for refining intervals or disks, that are already known 
to be isolating for the roots of a polynomial. Eor large L, that is if L dominates 
other parameters, their bit complexity is 0{nL). In comparison, using CIsolate 
for the refinement directly, its bit complexity would be of size O^n^L). We 
suspect that this bound can be further improved by using a proper modification 
of the Tl.-test, which only needs to evaluate F and its first derivative, and 
approximate multipoint evaluation. We have not analyzed these extensions, 
however, we are confident that our approach yields similar bit complexity bounds 
as provided in ]28] for the modified variant of Pan’s method. This will be subject 
of future work. 

Personal communication with Alexander Kobel from Max-Planck-Institute for Informatics 
in Saarbriicken. 
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7 Appendix 

7.1 Missing proofs in Section 3.1 

We split the proof of Theorem 2 into two technical lemmas: 

Lemma 15. Let A := A(m,r) be a disk that is (l,4c2 • maxi(fc) • n^)-isolating 
for the roots zi,... then, for all z G C 2 n^ ■ A, it holds that 0. 

Furthermore, 


E 


F^'^\m){cin ■ rf ^kl 


F^^\m)i\ 


< 


1 

2K' 


Proof. 1. For the first part, we may assume that fc > 1. Then, for the fc’th 
derivative of F and any complex z that is not a root of F, it holds that 


i?(fe)(z) 
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By way of contradiction, assume F^^\z) = 0 for some z £ C2U^ • A. Then, 
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Assuming k < n/2, we proceed to show 
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< Y,{i/2)>^-'^'/{k-ky. 

fc'=0 

< — 1 < 1, a contradiction. 

The preceding argument assumes k < n/2 because this allows us to freely 
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choose J and J' as indicated. But suppose k > nj^. We then have 
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2. Similar as above, with ..., denoting the roots of it holds 
that 
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where we used (5) for the last inequality. □ 

Lemma 16. Le X he a real value with X > 16K ■ maxi(fc)^ • n and suppose that 
A := A(m,r) is a disk that is {1, X)-isolating for the roots zi,...,Zk of F, then 


E 

i<k 


(m)| (cin • r)® ^k\ 
|F('®)(m)| i\ 
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1 


Proof. We may assume that fc > 1. Write F{x) = G{x)H{x) with G{x) = 
nLi(2^ - Zi) and H{x) = nj=/c+i {x — Zj). By induction, one shows that 
= EU and = fc!E,e(H) = 

k\ ■ X]jg( ["I ) ~ Zj)- It follows that 
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For G, we have G^''\x) = *! Ilj^ j {x — Zj), and thus |G'*^*^(to)| 
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Hence, using that A > l&Kk^ ■ n and cin > , ^ i \ , yields 
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□ 


7.2 Missing Proofs in Section 3.2 


Theorem (Restatement of Theorem 3). Denote the roots of F by zi,... ,Zn, 
then it holds that F[^l(a:) = '}2'i=o • nr=i(^ ~ ^i)- particular, the 

roots of the first Graeffe iterate are the squares of the roots of F. In addition, 
we have 


• maxidlFlUf > ||fW||^ > ||F||L • 2-4". 

Proof. Notice that On^ = follows directly from the definition of Further¬ 
more, we have 


FW(zf) = (-1)" . [F^izff - zf ■ F^izff] 

= {-ir ■ [Fe{z!) - z. ■ F,{z^)] ■ [F,{zf) + z, ■ F„{zf)] 
= (-!)" •[F,(z2)-Z4-F„(z2)].F(z4)=0. 


Going from F to an arbitrary small perturbation F (which has only simple roots 
and for which zf y z| for all pairs of distinct roots Zi and z^ of F\ we conclude 
that each root of has multiplicity mult(z?, mult(z,-,F). 

Hence, the first claim follows. For the second claim, notice that the left inequality 
follows immediately from the fact that each coefficient of is the sum of at 
most rif many products of the form zboi • aj, and each of these products has 
absolute value smaller than or equal to maxi(||F'||oo)^- For the right inequality 
we have to work harder: W.l.o.g., we may assume that \zi\ <2 for i = 1,..., fc, 
and that \zi\ > 2 for i = k + 1,... ,n. Let ^max be a point in the closure of 
the unit disk A(0,1) such that |F(zinax)| = maxa.| 2 |<i \F{z)\. Since F takes 
its maximum on the boundary of A(0,1), we must have [zmaxl = !> and using 
Cauchy’s Integral Theorem to write the coefficients of F in terms of an integral. 
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we conclude that |-F(zmax)l > Halloo- In addition, it holds that |-F(zmax)| < 
Yh=o |n*l • kmaxl” = Yh=o |a»| < (n + 1) • || J^||oo, and thus 

lll^lloo < |i"(w)| = max \F{z)\ < (n + 1) • \\F\\^. (27) 

z:\z\<l 


Applying the latter result to the polynomial g{x) := 11^=1 (-^ ~ ^i) yields the 
existence of a point z' with \z'\ = 1 and |( 7 (- 2 :')| > 1. Hence, it follows that 


|fW(z')| = |a„p |z'- \z'- zf\> \an\^ ]J z'- Zi) ■ {Vf + Zi)\ 
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where we used that |a; — j/| <3 for arbitrary complex points x,y with \x\ = 1 
and |j/| < 2, and that |a; — 2 | > for arbitrary complex points x,y,z with 

kl = ll/l = 1 and \z\ > 2. We conclude that 
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